///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/rest_gem_pow.cpp
///
/// \brief   implementation file for the restricted geminal power coefficient object
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<vector>
#include<complex>
#include<string>
#include<algorithm>
#include<cmath>

#include<formic/fqmc/rest_gem_pow.h>
#include<formic/fqmc/all_pair_jast.h>
#include<formic/wfn/configuration.h>
#include<formic/fqmc/rgp_detail.h>
#include<formic/exception.h>
#include<formic/all_data.h>
#include<formic/command/command.h>
#include<formic/mpi/interface.h>
#include<formic/lapack/interface.h>
#include<formic/numeric/numeric.h>
#include<formic/numeric/matrix.h>
#include<formic/output/output.h>
#include<formic/timing/timing.h>
#include<formic/vector/vector.h>
#include<formic/vector/matrix.h>
#include<formic/ham/ham.h>
#include<formic/ham/tei.h>

//#include<formic/fqmc/rgp_fortran.h>
#include<formic/fqmc/fortran/rgp_energy.h>
#include<formic/fqmc/fortran/rgp_der_hc.h>
#include<formic/fqmc/fortran/rgp_loc_en_adj_der.h>
#include<formic/fqmc/fortran/jrgp_loc_en_adj_der.h>
#include<formic/fqmc/fortran/jrgp_loc_en_integral_der.h>

namespace formic { namespace fqmc {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  A class for a restricted geminal power wavefunction
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class RGemPow : public formic::CoeffObj<S> {

    // private types
    private:

      /// \brief short name for a matrix type
      typedef formic::ConstMatrix<S,int> ConstMat;
      typedef formic::Matrix<S,int> Mat;

    // private data members
    private:

      /// \brief number of orbitals
      int _norb;

      /// \brief number of alpha electrons
      int _nap;

      /// \brief pointer to the pairing matrix
      S * _pm;

      /// \brief whether the pairing matrix is fixed or variable
      bool _fixed_pm;

      /// \brief whether to compute the der_hc function via the adjoint code
      bool _der_hc_by_adj;

      /// \brief vector to hold pairing matrix data (not always used)
      std::vector<S> _pm_data;

      /// \brief the occupied pairing matrix
      std::vector<S> _opm;

      /// \brief the inverse of the occupied pairing matrix
      std::vector<S> _oim;

      /// \brief used to hold the contraction of the inverse matrix with derivative coefficients
      std::vector<S> _oic;

      /// \brief workspace
      std::vector<S> _swork;

      /// \brief integer workspace
      std::vector<int> _iwork;

      /// \brief a reference to the Hamiltonian object
      formic::Hamiltonian<S> & _ham;

      /// \brief whether the object is initialized
      bool _initialized;

      /// \brief the energy is penalized by _s_sqr_penalty * S^2
      double _s_sqr_penalty;

      /// \brief the configuration for which the object is initialized
      formic::Configuration _config;

      /// \brief the current configuration's coefficient 
      formic::MantExp<S> _curr_coeff;

      /// \brief pointer to an object that manages the all pair jastrow factors
      formic::fqmc::AllPairJast<S> * _apj;

      /// \brief shared pointer holding all pair jastrows (not always used)
      boost::shared_ptr<formic::fqmc::AllPairJast<S> > _apj_shared_ptr;

      // tensors used for local energy and hamiltonian contraction
      std::vector<S> _Tia_aa;       ///<  alpha one electron excitation integrals
      std::vector<S> _Tia_bb;       ///<  beta  one electron excitation integrals
      std::vector<S> _Tiajb_aaaa;   ///<  all-alpha  two electron excitation integrals
      std::vector<S> _Tiajb_bbbb;   ///<  all-beta   two electron excitation integrals
      std::vector<S> _Tiajb_aabb;   ///<  alpha-beta two electron excitation integrals
      std::vector<S> _Bkl;          ///<  inverse matrix
      std::vector<S> _Blk;          ///<  transpose of inverse matrix
      std::vector<S> _Via;          ///<  row matrix for alpha excitations, second index labels row
      std::vector<S> _Ujb;          ///<  col matrix for beta  excitations, second index labels col
      std::vector<S> _Pab;          ///<  unoccupied pairing matrix
      std::vector<S> _VBia;         ///<  contraction of V and B
      std::vector<S> _BUjb;         ///<  contraction of B and U
      std::vector<S> _TVBjb_aaaa;   ///<  contraction of T, V, and B
      std::vector<S> _TVBjb_bbbb;   ///<  contraction of T, U, and B
      std::vector<S> _TVBjb_aabb;   ///<  contraction of T, V, and B
      std::vector<S> _TVBaj;        ///<  contraction of T, V, and B
      std::vector<S> _TUBaj;        ///<  contraction of T, U, and B
      std::vector<S> _VBUab;        ///<  contraction of V, B, and U
      std::vector<S> _TBab;         ///<  contraction of T and B

      // more tensors used for hamiltonian contraction
      std::vector<S> _Ckl;          ///<  derivative contracted inverse matrix
      std::vector<S> _Clk;          ///<  transpose of derivative contracted inverse matrix
      std::vector<S> _Zia;          ///<  row matrix for alpha excitations built from derivatives
      std::vector<S> _Yjb;          ///<  col matrix for beta  excitations built from derivatives
      std::vector<S> _Qab;          ///<  unoccupied pairing matrix built from derivatives
      std::vector<S> _ZBia;         ///<  contraction of Z and B
      std::vector<S> _VCia;         ///<  contraction of V and C
      std::vector<S> _BYjb;         ///<  contraction of B and Y
      std::vector<S> _CUjb;         ///<  contraction of C and U
      std::vector<S> _dVBUab;       ///<  derivative contracted VBU
      std::vector<S> _TCab;         ///<  contraction of T and C
      std::vector<S> _TBUia;        ///<  contraction of T, B, and U

      // tensors used for adjoint calculation
      std::vector<S> _RTA;          ///<  adjoint for unoccupied rows times inverse matrix
      std::vector<S> _Ksi;          ///<  adjoint for single index occ jastrow intermediate
      std::vector<S> _Ksj;          ///<  adjoint for single index occ jastrow intermediate
      std::vector<S> _Ksa;          ///<  adjoint for single index unocc jastrow intermediate
      std::vector<S> _Ksb;          ///<  adjoint for single index unocc jastrow intermediate
      std::vector<S> _Kdia;         ///<  adjoint for double index occ-unocc jastrow intermediate
      std::vector<S> _Kdij;         ///<  adjoint for double index occ-occ jastrow intermediate
      std::vector<S> _Kdib;         ///<  adjoint for double index occ-unocc jastrow intermediate
      std::vector<S> _Kdaj;         ///<  adjoint for double index unocc-occ jastrow intermediate
      std::vector<S> _Kdab;         ///<  adjoint for double index unocc-unocc jastrow intermediate
      std::vector<S> _Kdjb;         ///<  adjoint for double index occ-unocc jastrow intermediate

    // private member functions
    private:

      /// \brief disable the assignment operator
      RGemPow<S> & operator=(const RGemPow<S> & other);

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  copy constructor is private and intended to be used by the clone function
      ///
      /// \param[in]      other    the object to copy from
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      RGemPow(const RGemPow<S> & other)
        : formic::CoeffObj<S>(other),
          _norb(other._norb),
          _nap(other._nap),
          _pm(other._pm),
          _fixed_pm(other._fixed_pm),
          _der_hc_by_adj(other._der_hc_by_adj),
          _pm_data(other._pm_data),
          _opm(other._opm),
          _oim(other._oim),
          _oic(other._oic),
          _swork(other._swork),
          _iwork(other._iwork),
          _ham(other._ham),
          _initialized(other._initialized),
          _s_sqr_penalty(other._s_sqr_penalty),
          _config(other._config),
          _curr_coeff(other._curr_coeff),
          _apj(other._apj),
          _apj_shared_ptr(other._apj_shared_ptr),
          _Tia_aa(other._Tia_aa),
          _Tia_bb(other._Tia_bb),
          _Tiajb_aaaa(other._Tiajb_aaaa),
          _Tiajb_bbbb(other._Tiajb_bbbb),
          _Tiajb_aabb(other._Tiajb_aabb),
          _Bkl(other._Bkl),
          _Blk(other._Blk),
          _Via(other._Via),
          _Ujb(other._Ujb),
          _Pab(other._Pab),
          _VBia(other._VBia),
          _BUjb(other._BUjb),
          _TVBjb_aaaa(other._TVBjb_aaaa),
          _TVBjb_bbbb(other._TVBjb_bbbb),
          _TVBjb_aabb(other._TVBjb_aabb),
          _TVBaj(other._TVBaj),
          _TUBaj(other._TUBaj),
          _VBUab(other._VBUab),
          _TBab(other._TBab),
          _Ckl(other._Ckl),
          _Clk(other._Clk),
          _Zia(other._Zia),
          _Yjb(other._Yjb),
          _Qab(other._Qab),
          _ZBia(other._ZBia),
          _VCia(other._VCia),
          _BYjb(other._BYjb),
          _CUjb(other._CUjb),
          _dVBUab(other._dVBUab),
          _TCab(other._TCab),
          _TBUia(other._TBUia),
          _RTA(other._RTA),
          _Ksi(other._Ksi),
          _Ksj(other._Ksj),
          _Ksa(other._Ksa),
          _Ksb(other._Ksb),
          _Kdia(other._Kdia),
          _Kdij(other._Kdij),
          _Kdib(other._Kdib),
          _Kdaj(other._Kdaj),
          _Kdab(other._Kdab),
          _Kdjb(other._Kdjb)
      {

        // prepare the pairing matrix (necessary if other's _pm points to the wavefunction)
        _pm_data.resize(_norb*_norb);
        std::copy(other._pm, other._pm + _norb*_norb, _pm_data.begin());
        _pm = &_pm_data.at(0);

        // clone the all pair jastrow object
        typedef formic::fqmc::AllPairJast<S> typeA;
        typedef formic::CoeffObj<S> typeC;
        _apj_shared_ptr.reset();
        if (other._apj)
          _apj_shared_ptr = boost::dynamic_pointer_cast<typeA, typeC>(other._apj->clone());
        _apj = _apj_shared_ptr.get();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief checks that the object has been initialized for a configuration
      ///
      /// \param[in]      calling_func  the name of the function in which this check is performed
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void check_initialized(const std::string & calling_func) {

        if ( !_initialized )
          throw formic::Exception("%s requires the object to first be initialized") % calling_func;

      }

      /// \brief  position of the pairing matrix in the derivatives ordering
      int pm_pos(const bool ign_fixed = true) const { return 0; }

      /// \brief  position of the all pair jastrow factor variables in the derivatives ordering
      int jf_pos(const bool ign_fixed = true) const { return ( ign_fixed && _fixed_pm ? 0 : _norb*_norb ); }

//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief checks that the supplied configuration difference object is compatible with the stored
//      ///        configuration
//      ///
//      /// \param[in]      diff          object describing a change to the configuration
//      /// \param[in]      calling_func  the name of the function in which this check is performed
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void check_config_diff(const formic::ConfigDiff & diff, const std::string & calling_func) {
//
//        // make sure the object is initialized
//        this->check_initialized("RGemPow::check_config_diff");
//
//        // check newly full sites
//        for (int i = 0; i < diff.newly_full().size(); i++)
//          if ( _config[diff.newly_full().at(i)] != 0 )
//            throw formic::Exception("cannot move an electron into a site that is already occupied in %s") % calling_func;
//
//        // check newly empty sites
//        for (int i = 0; i < diff.newly_empty().size(); i++)
//          if ( _config[diff.newly_empty().at(i)] == 0 )
//            throw formic::Exception("cannot move an electron out of a site that is empty in %s") % calling_func;
//
//      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Prepares the tensors of integrals necessary for computing the alpha-alpha and
      ///        alpha-alpha-alpha-alpha contributions to the local energy.
      ///
      /// \param[in]      ao2s          alpha map from occupied to site
      /// \param[in]      bo2s          beta  map from occupied to site
      /// \param[in]      au2s          alpha map from unoccupied to site
      /// \param[in]      bu2s          beta  map from unoccupied to site
      /// \param[out]     Tia           vector to hold the resulting one electron integrals
      /// \param[out]     Tiajb         vector to hold the resulting two electron integrals
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void prep_local_energy_ints_aa_aaaa(const std::vector<int> & ao2s,
                                          const std::vector<int> & bo2s,
                                          const std::vector<int> & au2s,
                                          const std::vector<int> & bu2s,
                                          std::vector<S> & Tia,
                                          std::vector<S> & Tiajb)
      {

        // get number of occupied and unoccupied orbitals
        const int no = _nap;
        const int nv = _norb - _nap;

        // get pointer to dense two electron integrals
        const S * const all_tei = _ham.dense_tei_ptr();

        // size the integral tensors
        if ( Tia.size() != no*nv )
          Tia.resize(no*nv);
        if ( Tiajb.size() != no*nv*no*nv )
          Tiajb.resize(no*nv*no*nv);

        // prepare one electron integral tensor Tia

        // ia one electron integral
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
          Tia[i+no*a] = _ham.get_mei(ao2s[i], au2s[a]);

        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++) {

          const S integral = all_tei[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))];

          // (ii|ia) and (ii|ja), although the code exchanges i and j
          Tia[i+no*a] += 0.5 * integral;

          // (ia|jj)
          Tia[i+no*a] += 0.5 * integral;

        }

        // (ib|ba) and (ia|aa)
        for (int a = 0; a < nv; a++)
        for (int b = 0; b < nv; b++)
        for (int i = 0; i < no; i++)
          Tia[i+no*a] += 0.5 * all_tei[ao2s[i]+_norb*(au2s[b]+_norb*(au2s[b]+_norb*au2s[a]))];

        // (ia|ji), although the code exchanges i and j
        for (int a = 0; a < nv; a++)
        for (int j = 0; j < no; j++)
        for (int i = 0; i < no; i++)
          Tia[i+no*a] -= 0.5 * all_tei[ao2s[i]+_norb*(ao2s[j]+_norb*(ao2s[j]+_norb*au2s[a]))];

        // bbaa (jj|ia) and aabb (ia|jj)
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++)
          Tia[i+no*a] += all_tei[bo2s[j]+_norb*(bo2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))];

        // prepare two electron integrals
        formic::tei_get_ovov_aaaa_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), all_tei, &Tiajb.at(0));

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Prepares the tensors of integrals necessary for computing the alpha-alpha and
      ///        alpha-alpha-alpha-alpha contributions to the local energy.
      ///
      /// \param[in]      weight    weight to apply to the derivatives
      /// \param[in]      ao2s      alpha map from occupied to site
      /// \param[in]      bo2s      beta  map from occupied to site
      /// \param[in]      au2s      alpha map from unoccupied to site
      /// \param[in]      bu2s      beta  map from unoccupied to site
      /// \param[in]      Wia       derivative adjoints of occ-unocc one electron integrals
      /// \param[in]      Wiajb     derivative adjoints of occ-unocc-occ-unocc two electron integrals
      /// \param[in,out]  mei_adj   cumulative derivatives of all modified one electron integrals
      /// \param[in,out]  tei_adj   cumulative derivatives of all two electron integrals
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void add_local_energy_int_der_aa_aaaa(const S weight,
                                            const std::vector<int> & ao2s,
                                            const std::vector<int> & bo2s,
                                            const std::vector<int> & au2s,
                                            const std::vector<int> & bu2s,
                                            const std::vector<S> & Wia,
                                            const std::vector<S> & Wiajb,
                                            S * const mei_adj,
                                            S * const tei_adj)
      {

        // get number of occupied and unoccupied orbitals
        const int no = _nap;
        const int nv = _norb - _nap;

//        // get pointer to dense two electron integrals
//        const S * const all_tei = _ham.dense_tei_ptr();
//
//        // size the integral tensors
//        if ( Tia.size() != no*nv )
//          Tia.resize(no*nv);
//        if ( Tiajb.size() != no*nv*no*nv )
//          Tiajb.resize(no*nv*no*nv);

        // prepare one electron integral tensor Tia

        // ia one electron integral
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
          mei_adj[ ao2s[i] * _norb + au2s[a] ] += weight * Wia[i+no*a];
          //Tia[i+no*a] = _ham.get_mei(ao2s[i], au2s[a]);
          //T get_mei(const int p, const int q) const { return _mei[ p * _norb + q ]; }

        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++) {

          S & integral = tei_adj[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))];
          //const S integral = all_tei[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))];

          // (ii|ia) and (ii|ja), although the code exchanges i and j
          integral += 0.5 * weight * Wia[i+no*a];
          //Tia[i+no*a] += 0.5 * integral;

          // (ia|jj)
          integral += 0.5 * weight * Wia[i+no*a];
          //Tia[i+no*a] += 0.5 * integral;

        }

        // (ib|ba) and (ia|aa)
        for (int a = 0; a < nv; a++)
        for (int b = 0; b < nv; b++)
        for (int i = 0; i < no; i++)
          tei_adj[ao2s[i]+_norb*(au2s[b]+_norb*(au2s[b]+_norb*au2s[a]))] += weight * 0.5 * Wia[i+no*a];
          //Tia[i+no*a] += 0.5 * all_tei[ao2s[i]+_norb*(au2s[b]+_norb*(au2s[b]+_norb*au2s[a]))];

        // (ia|ji), although the code exchanges i and j
        for (int a = 0; a < nv; a++)
        for (int j = 0; j < no; j++)
        for (int i = 0; i < no; i++)
          tei_adj[ao2s[i]+_norb*(ao2s[j]+_norb*(ao2s[j]+_norb*au2s[a]))] -= weight * 0.5 * Wia[i+no*a];
          //Tia[i+no*a] -= 0.5 * all_tei[ao2s[i]+_norb*(ao2s[j]+_norb*(ao2s[j]+_norb*au2s[a]))];

        // bbaa (jj|ia) and aabb (ia|jj)
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++)
          tei_adj[bo2s[j]+_norb*(bo2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))] += weight * Wia[i+no*a];
          //Tia[i+no*a] += all_tei[bo2s[j]+_norb*(bo2s[j]+_norb*(ao2s[i]+_norb*au2s[a]))];

        // prepare two electron integrals
        const int nt = no + nv;
        for (int b = 0; b < nv; b++)
        for (int j = 0; j < no; j++)
        for (int a = 0; a < b  ; a++)
        for (int i = 0; i < j  ; i++) {
          S & integ1 = tei_adj[ ao2s[i] + nt * ( au2s[a] + nt * ( ao2s[j] + nt * au2s[b] ) ) ];
          S & integ2 = tei_adj[ ao2s[j] + nt * ( au2s[a] + nt * ( ao2s[i] + nt * au2s[b] ) ) ];
          integ1 += weight * Wiajb[ i + no * ( a + nv * ( j + no * b ) ) ];
          integ1 += weight * Wiajb[ j + no * ( b + nv * ( i + no * a ) ) ];
          integ2 += weight * Wiajb[ j + no * ( a + nv * ( i + no * b ) ) ];
          integ2 += weight * Wiajb[ i + no * ( b + nv * ( j + no * a ) ) ];
        }
        //const int nt = no + nv;
        //std::fill(tei_slice, tei_slice + no*no*nv*nv, formic::zero(S()));
        //for (int b = 0; b < nv; b++)
        //for (int j = 0; j < no; j++)
        //for (int a = 0; a < b  ; a++)
        //for (int i = 0; i < j  ; i++) {
        //  const S integ1 = tei_dense[ ao2s[i] + nt * ( au2s[a] + nt * ( ao2s[j] + nt * au2s[b] ) ) ];
        //  const S integ2 = tei_dense[ ao2s[j] + nt * ( au2s[a] + nt * ( ao2s[i] + nt * au2s[b] ) ) ];
        //  tei_slice[ i + no * ( a + nv * ( j + no * b ) ) ] = integ1;
        //  tei_slice[ j + no * ( b + nv * ( i + no * a ) ) ] = integ1;
        //  tei_slice[ j + no * ( a + nv * ( i + no * b ) ) ] = integ2;
        //  tei_slice[ i + no * ( b + nv * ( j + no * a ) ) ] = integ2;
        //}
        //formic::tei_get_ovov_aaaa_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), all_tei, &Tiajb.at(0));

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Prepares the new row, new column, and unoccupied pairing matrices for the current
      ///        configuration using the supplied full pairing matrix
      ///
      /// \param[in]      pm            the full pairing matrix
      /// \param[out]     Via           on exit, the new row matrix
      /// \param[out]     Ujb           on exit, the new column matrix
      /// \param[out]     Pab           on exit, the unoccupied pairing matrix
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void prep_row_col_unocc_mat(const S * const pm,
                                  std::vector<S> & Via,
                                  std::vector<S> & Ujb,
                                  std::vector<S> & Pab)
      {

        // get number of occupied and unoccupied orbitals
        const int no = _nap;
        const int nv = _norb - _nap;

        // get references to occupation arrays
        const std::vector<int> & ao2s = _config.a_occ_to_site();
        const std::vector<int> & bo2s = _config.b_occ_to_site();
        const std::vector<int> & au2s = _config.a_unocc_to_site();
        const std::vector<int> & bu2s = _config.b_unocc_to_site();

        // prepare new row matrix
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
          Via[i+no*a] = pm[ au2s[a] * _norb + bo2s[i] ];

        // prepare new column matrix
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++)
          Ujb[i+no*a] = pm[ ao2s[i] * _norb + bu2s[a] ];

        // prepare unoccupied pairing matrix
        for (int a = 0; a < nv; a++)
        for (int b = 0; b < nv; b++)
          Pab[a+nv*b] = pm[ au2s[a] * _norb + bu2s[b] ];

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Returns the non-excitation contributions to the local energy from Hamiltonian terms
      ///         with aa, aaaa, and aabb integrals.
      ///
      /// \param[in]      ao2s          alpha map from occupied to site
      /// \param[in]      bo2s          beta  map from occupied to site
      /// \param[in]      au2s          alpha map from unoccupied to site
      /// \param[in]      bu2s          beta  map from unoccupied to site
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S non_excite_energy_aa_aaaa_aabb(const std::vector<int> & ao2s,
                                       const std::vector<int> & bo2s,
                                       const std::vector<int> & au2s,
                                       const std::vector<int> & bu2s)
      {

        // get number of occupied and unoccupied orbitals
        const int no = _nap;
        const int nv = _norb - _nap;

        // get pointer to dense two electron integrals
        const S * const all_tei = _ham.dense_tei_ptr();

        // initialize non-excitation contribution
        S retval = formic::zero(S());

        // aa
        for (int i = 0; i < no; i++)
          retval += _ham.get_mei(ao2s[i], ao2s[i]);

        // aaaa (ii|ii) and (ii|jj)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++)
          retval += 0.5 * all_tei[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*ao2s[i]))];

        // aaaa (ia|ai)
        for (int i = 0; i < no; i++)
        for (int a = 0; a < nv; a++)
          retval += 0.5 * all_tei[au2s[a]+_norb*(ao2s[i]+_norb*(au2s[a]+_norb*ao2s[i]))];

        // aabb (ii|jj) non-excitations
        for (int j = 0; j < no; j++)
        for (int i = 0; i < no; i++)
          retval += 0.5 * all_tei[ao2s[i]+_norb*(ao2s[i]+_norb*(bo2s[j]+_norb*bo2s[j]))];

        // penalty for non-zero s_squared
        if ( _s_sqr_penalty != 0.0 ) {
          retval += formic::unity(S()) * _s_sqr_penalty * ( 0.5 * 0.25 * ( ao2s.size() - bo2s.size() ) * ( ao2s.size() - bo2s.size() ) );
          assert( 2*_norb == _config.fock_config().size() );
          for (int i = 0; i < _norb; i++)
            if ( _config[2*i+0] != _config[2*i+1] )
              retval += _s_sqr_penalty * 0.5 * 0.5 * formic::unity(S());
        }

        // return the non excitation contribution
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Adds the non-excitation related local energy integral derivative contributions.
      ///
      /// \param[in]      ao2s      alpha map from occupied to site
      /// \param[in]      bo2s      beta  map from occupied to site
      /// \param[in]      au2s      alpha map from unoccupied to site
      /// \param[in]      bu2s      beta  map from unoccupied to site
      /// \param[in]      weight    weight to apply to the derivatives
      /// \param[in,out]  mei_adj   cumulative adjoint derivatives for modiefied one electron integrals
      /// \param[in,out]  tei_adj   cumulative adjoint derivatives for two electron integrals
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void non_excite_en_integral_der_aa_aaaa_aabb(const std::vector<int> & ao2s,
                                                   const std::vector<int> & bo2s,
                                                   const std::vector<int> & au2s,
                                                   const std::vector<int> & bu2s,
                                                   const S weight,
                                                   S * const mei_adj,
                                                   S * const tei_adj)
      {

        // get number of occupied and unoccupied orbitals
        const int no = _nap;
        const int nv = _norb - _nap;

        //// get pointer to dense two electron integrals
        //const S * const all_tei = _ham.dense_tei_ptr();

        //// initialize non-excitation contribution
        //S retval = formic::zero(S());

        // aa
        for (int i = 0; i < no; i++)
          mei_adj[ ao2s[i] * _norb + ao2s[i] ] += weight;
          //retval += _ham.get_mei(ao2s[i], ao2s[i]);
          //T get_mei(const int p, const int q) const { return _mei[ p * _norb + q ]; }

        // aaaa (ii|ii) and (ii|jj)
        for (int i = 0; i < no; i++)
        for (int j = 0; j < no; j++)
          tei_adj[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*ao2s[i]))] += 0.5 * weight;
          //retval += 0.5 * all_tei[ao2s[j]+_norb*(ao2s[j]+_norb*(ao2s[i]+_norb*ao2s[i]))];

        // aaaa (ia|ai)
        for (int i = 0; i < no; i++)
        for (int a = 0; a < nv; a++)
          tei_adj[au2s[a]+_norb*(ao2s[i]+_norb*(au2s[a]+_norb*ao2s[i]))] += 0.5 * weight;
          //retval += 0.5 * all_tei[au2s[a]+_norb*(ao2s[i]+_norb*(au2s[a]+_norb*ao2s[i]))];

        // aabb (ii|jj) non-excitations
        for (int j = 0; j < no; j++)
        for (int i = 0; i < no; i++)
          tei_adj[ao2s[i]+_norb*(ao2s[i]+_norb*(bo2s[j]+_norb*bo2s[j]))] += 0.5 * weight;
          //retval += 0.5 * all_tei[ao2s[i]+_norb*(ao2s[i]+_norb*(bo2s[j]+_norb*bo2s[j]))];

        //// return the non excitation contribution
        //return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  computes the derivative rations for the current configuration, optionally skipping the
      ///         jastrow ratios
      ///
      /// \param[in]      do_jr      whether to do the jastrow ratios
      /// \param[out]     ratios     on exit, the derivative ratios.
      ///                            for ign_fixed == true, size == this->n_der()
      ///                            for ign_fixed == false, size == this->n_var()
      /// \param[in]      ign_fixed  whether to ignore variables that are fixed
      ///
      ///   The ratio for variable i is     <n|Psi_i> / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void der_ratios_detail(const bool do_jr, S * const ratios, const bool ign_fixed = true) {

        // get references to occupation arrays
        const std::vector<int> & ao2s = _config.a_occ_to_site();
        const std::vector<int> & bo2s = _config.b_occ_to_site();

        // do pairing matrix ratios
        if ( !_fixed_pm || !ign_fixed ) {
          S * const pm_ratios = ratios + this->pm_pos(ign_fixed);
          std::fill(pm_ratios, pm_ratios + _norb*_norb, formic::zero(S()));
          for (int i = 0; i < _nap; i++) {
            const int ai = ao2s[i];
            for (int j = 0; j < _nap; j++) {
              const int bj = bo2s[j];
              assert( _norb * ai + bj <  _norb*_norb );
              assert( _norb * ai + bj >= 0           );
              assert( _nap * i + j <  _oim.size() );
              assert( _nap * i + j >= 0           );
              pm_ratios[ _norb * ai + bj ] = _oim[ _nap * j + i ];
            }
          }
        }

        // do jastrow ratios
        if ( _apj && do_jr )
          _apj->der_ratios(ratios + this->jf_pos(ign_fixed), ign_fixed);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Adds the unoccupied row matrix's adjoint derivative contribution to the
      ///         supplied vector of local energy derivatives.
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      au2s     vector of unoccupied alpha sites
      /// \param[in]      bo2s     vector of occupied beta sites
      /// \param[in]      urm_adj  adjoint derivatives of local energy w.r.t. the unoccupied row matrix
      /// \param[in,out]  der_out  the vector of local energy derivatives w.r.t. pairing matrix variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void urm_adjoint_to_der_vec(const int no,
                                  const int nv,
                                  const std::vector<int> & au2s,
                                  const std::vector<int> & bo2s,
                                  const ConstMat & urm_adj, // (nv,no)  (alpha, beta)
                                  S * const der_out) {

        // expand unoccupied row derivatives
        for (int j = 0; j < no; j++) // beta
        for (int a = 0; a < nv; a++) // alpha
          der_out[ this->pm_pos() + au2s[a] * (no+nv) + bo2s[j] ] += urm_adj(a,j);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Adds the inverse-occupied-pairing-matrix's adjoint derivative contribution to the
      ///         supplied vector of local energy derivatives.
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      ao2s     vector of occupied alpha sites
      /// \param[in]      bo2s     vector of occupied beta sites
      /// \param[in]      oim_adj  adjoint derivatives of local energy w.r.t. the inverse-occupied-matrix
      /// \param[in]      cit      transpose of the inverse-occupied-pairing-matrix
      /// \param[in,out]  der_out  the vector of local energy derivatives w.r.t. pairing matrix variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void ipm_adjoint_to_der_vec(const int no,
                                  const int nv,
                                  const std::vector<int> & ao2s,
                                  const std::vector<int> & bo2s,
                                  const ConstMat & oim_adj, // (no,no)  (beta, alpha)
                                  const ConstMat & cit, // transpose of inverse pairing matrix
                                  S * const der_out) {

        // get adjoint for occupied pairing matrix
        const Mat opm_adj = -( cit * oim_adj ) * cit;

        // expand occupied pairing matrix derivatives
        for (int j = 0; j < no; j++) // beta
        for (int i = 0; i < no; i++) // alpha
          der_out[ this->pm_pos() + ao2s[i] * (no+nv) + bo2s[j] ] += opm_adj(i,j);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Adds the RT matrix's adjoint derivative contribution to the
      ///         supplied vector of local energy derivatives.
      ///         (The RT matrix is the product of the unoccupied row matrix with the
      ///          inverse-occupied-pairing-matrix.)
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      ao2s     vector of occupied alpha sites
      /// \param[in]      bo2s     vector of occupied beta sites
      /// \param[in]      au2s     vector of unoccupied alpha sites
      /// \param[in]      bu2s     vector of unoccupied beta sites
      /// \param[in]      rta      adjoint derivatives of local energy w.r.t. the RT matrix
      /// \param[in]      cur_im   the current inverse of the occupied pairing matrix
      /// \param[in]      cur_urm  the current unoccupied row matrix
      /// \param[in,out]  der_out  the vector of local energy derivatives w.r.t. pairing matrix variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void rtheta_adjoint_to_der_vec(const int no,
                                     const int nv,
                                     const std::vector<int> & ao2s,
                                     const std::vector<int> & bo2s,
                                     const std::vector<int> & au2s,
                                     const std::vector<int> & bu2s,
                                     const ConstMat & rta,
                                     const ConstMat & cur_im,
                                     const ConstMat & cur_urm,
                                     S * const der_out) {

        // get transpose of current inverse matrix
        const Mat cit = cur_im.t();

        // get adjoint for unoccupied row matrix ( dimension nv by no )
        const Mat uradj = rta * cit;

        // expand unoccupied row derivatives
        this->urm_adjoint_to_der_vec(no, nv, au2s, bo2s, uradj, der_out);

        // get adjoint for inverse pairing matrix
        const Mat oim_adj = cur_urm.t() * rta;

        // get derivatives resulting from oim_adj as adjoint for inverse pairing matrix
        this->ipm_adjoint_to_der_vec(no, nv, ao2s, bo2s, oim_adj, cit, der_out);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Adds the unoccupied pairing matrix's adjoint derivative contribution to the
      ///         supplied vector of local energy derivatives.
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      au2s     vector of unoccupied alpha sites
      /// \param[in]      bu2s     vector of unoccupied beta sites
      /// \param[in]      upm_adj  adjoint derivatives of local energy w.r.t. the unoccupied pairing matrix
      /// \param[in,out]  der_out  the vector of local energy derivatives w.r.t. pairing matrix variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void upm_adjoint_to_der_vec(const int no,
                                  const int nv,
                                  const std::vector<int> & au2s,
                                  const std::vector<int> & bu2s,
                                  const ConstMat & upm_adj, // (nv,nv)  (alpha, beta)
                                  S * const der_out) {

        // expand unoccupied pairing matrix derivatives
        for (int b = 0; b < nv; b++) // beta
        for (int a = 0; a < nv; a++) // alpha
          der_out[ this->pm_pos() + au2s[a] * (no+nv) + bu2s[b] ] += upm_adj(a,b);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the aa one-electron integral and aaaa two-electron integral contributions to
      ///         the local energy and its derivatives w.r.t. pairing matrix and jastrow variables.
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      ao2s     vector of occupied alpha sites
      /// \param[in]      bo2s     vector of occupied beta sites
      /// \param[in]      au2s     vector of unoccupied alpha sites
      /// \param[in]      bu2s     vector of unoccupied beta sites
      /// \param[in]      cur_im   the current inverse of the occupied pairing matrix
      /// \param[in]      cur_pm   the current all-site pairing matrix
      /// \param[in]      Jsi      occ jastrow intermediate
      /// \param[in]      Jsa      unocc jastrow intermediate
      /// \param[in]      Jdia     occ-unocc jastrow intermediate
      /// \param[in]      Jdai     unocc-occ jastrow intermediate (just the transpose of the Jdia matrix)
      /// \param[in]      Jdij     occ-occ jastrow intermediate
      /// \param[in]      Jdab     unocc-unocc jastrow intermediate
      /// \param[in,out]  Tia      vector used to store the needed slice of the one electron integrals
      /// \param[in,out]  Tiajb    vector used to store the needed slice of the two electron integrals
      /// \param[in,out]  local_e  the local energy
      /// \param[in,out]  der_out  the vector of local energy derivatives
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void loc_en_der_aa_aaaa(const int no,
                              const int nv,
                              const std::vector<int> & ao2s,
                              const std::vector<int> & bo2s,
                              const std::vector<int> & au2s,
                              const std::vector<int> & bu2s,
                              const ConstMat & cur_im,
                              const ConstMat & cur_pm,
                              const S * const  Jsi,
                              const S * const  Jsa,
                              const S * const Jdia,
                              const S * const Jdai,
                              const S * const Jdij,
                              const S * const Jdab,
                              std::vector<S> & Tia,
                              std::vector<S> & Tiajb,
                              S & local_e,
                              S * const der_out) {

        // prepare bare aa and aaaa integrals
        this->prep_local_energy_ints_aa_aaaa(ao2s, bo2s, au2s, bu2s, Tia, Tiajb);

        // get the unoccupied row elements of the current pairing matrix
        const ConstMat current_urm = cur_pm.slice(au2s, bo2s);

        // prepare contraction of new row matrix with inverse of occupied pairing matrix
        const ConstMat rt = current_urm * cur_im;

        // get matrix to hold the derivative adjoint of the rt matrix (note indices are transposed for now)
        Mat rta(no,nv);

        // compute energy and adjoints
        S ad_energy = formic::zero(S());
        //formic::start_timer("der_hc_autodiff_fortran");
        if (_apj) {
          formic::fqmc::jagp_local_energy_adjoints_aa_aaaa(no, nv, &rt(0,0), &rt.t()(0,0),
                                                           Jsi, Jsa, Jdia, Jdai, Jdij, Jdab,
                                                           &Tia.at(0), &Tiajb.at(0), &ad_energy, &rta(0,0),
                                                           &_Ksi.at(0), &_Ksa.at(0), &_Kdia.at(0), &_Kdaj.at(0), &_Kdij.at(0), &_Kdab.at(0));
        } else {
          formic::fqmc::agp_local_energy_adjoints_aa_aaaa(no, nv, &rt(0,0), &rt.t()(0,0), &Tia.at(0), &Tiajb.at(0), &ad_energy, &rta(0,0));
        }
        //formic::stop_timer("der_hc_autodiff_fortran");

        // add energy contribution
        local_e += ad_energy;

        // un-transpose the rta matrix, giving it dimension (nv,no)
        rta.tip();

        // get derivatives resulting from rthetha
        this->rtheta_adjoint_to_der_vec(no, nv, ao2s, bo2s, au2s, bu2s, rta, cur_im, current_urm, der_out);

        // get derivatives resulting from jastrows
        if (_apj) {
          //formic::start_timer("der_hc_autodiff_unpack_jf");
          S * aa_pos;
          S * bb_pos;
          S * ab_pos;
          _apj->get_aa_bb_ab_position_pointers(der_out + this->jf_pos(), aa_pos, bb_pos, ab_pos);
          _apj->intermed_adj_to_der_vec_aa(&ao2s[0], &bo2s[0], &au2s[0], &bu2s[0],
                                           Jsi, Jsa, Jdij, Jdab, Jdia, Jdai,
                                           &_Ksi[0], &_Ksa[0], &_Kdij[0], &_Kdab[0], &_Kdia[0], &_Kdaj[0],
                                           aa_pos, ab_pos);
          //formic::stop_timer("der_hc_autodiff_unpack_jf");
        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the aabb two-electron integral contributions to the local energy and its
      ///         derivatives w.r.t. pairing matrix and jastrow variables.
      ///
      /// \param[in]      no         the number of occupied orbitals
      /// \param[in]      nv         the number of unoccupied orbitals
      /// \param[in]      ao2s       vector of occupied alpha sites
      /// \param[in]      bo2s       vector of occupied beta sites
      /// \param[in]      au2s       vector of unoccupied alpha sites
      /// \param[in]      bu2s       vector of unoccupied beta sites
      /// \param[in]      cur_im     the current inverse of the occupied pairing matrix
      /// \param[in]      cur_pm     the current all-site pairing matrix
      /// \param[in]      Jsi        jastrow single-index intermediate for index i
      /// \param[in]      Jsa        jastrow single-index intermediate for index a
      /// \param[in]      Jsj        jastrow single-index intermediate for index j
      /// \param[in]      Jsb        jastrow single-index intermediate for index b
      /// \param[in]      Jdia       jastrow double-index intermediate for indices i and a
      /// \param[in]      Jdij       jastrow double-index intermediate for indices i and j
      /// \param[in]      Jdib       jastrow double-index intermediate for indices i and b
      /// \param[in]      Jdaj       jastrow double-index intermediate for indices a and j
      /// \param[in]      Jdab       jastrow double-index intermediate for indices a and b
      /// \param[in]      Jdjb       jastrow double-index intermediate for indices j and b
      /// \param[in,out]  Tiajb      vector used to store the needed slice of the two electron integrals
      /// \param[in,out]  local_e    the local energy
      /// \param[in,out]  der_out_a  vector of local energy derivatives in alpha-first ordering
      /// \param[in,out]  der_out_b  vector of local energy derivatives in beta-first ordering
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void loc_en_der_aabb(const int no,
                           const int nv,
                           const std::vector<int> & ao2s,
                           const std::vector<int> & bo2s,
                           const std::vector<int> & au2s,
                           const std::vector<int> & bu2s,
                           const ConstMat & cur_im,
                           const ConstMat & cur_pm,
                           const S * const  Jsi,
                           const S * const  Jsa,
                           const S * const  Jsj,
                           const S * const  Jsb,
                           const S * const Jdia,
                           const S * const Jdij,
                           const S * const Jdib,
                           const S * const Jdaj,
                           const S * const Jdab,
                           const S * const Jdjb,
                           std::vector<S> & Tiajb,
                           S & local_e,
                           S * const der_out_a,
                           S * const der_out_b) {

        // prepare bare aabb integrals
        formic::tei_get_ovov_aabb_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), _ham.dense_tei_ptr(), &Tiajb.at(0));

        // get the unoccupied row elements of the current pairing matrix
        const ConstMat cur_urm = cur_pm.slice(au2s, bo2s);

        // get the unoccupied column elements of the current pairing matrix
        const ConstMat cur_ucm = cur_pm.slice(ao2s, bu2s);

        // get the unocc-unocc pairing matrix
        const ConstMat cur_uum = cur_pm.slice(au2s, bu2s);

        // prepare rt(a,i) = contraction of new row matrix with inverse of occupied pairing matrix
        const ConstMat rt = cur_urm * cur_im;

        // prepare tc(j,b) = contraction of new column matrix with inverse of occupied pairing matrix
        const ConstMat tc = cur_im * cur_ucm;

        // prepare rtc(a,b) = contraction of new row, inverse, and new column matrices
        const ConstMat rtc = cur_urm * tc;

        // need to get out VTC(a,i), VRT(j,b), VARTC(j,i), and VT(a,b)

        // get matrix to hold the VRT contraction
        Mat vrt(no,nv); // (beta,beta)

        // get matrix to hold the VTC contraction (note indices are transposed for now)
        Mat vtc(no,nv); // (alpha,alpha)

        // get matrix to hold the VT contraction
        Mat vt(nv,nv); // (alpha,beta)

        // get matrix to hold the VARTC contraction (note indices are transposed for now)
        Mat vartc(no,no); // (alpha,beta)

        // compute energy and adjoints
        S ad_energy = formic::zero(S());
        //formic::start_timer("der_hc_autodiff_fortran");
        if (_apj) {
          formic::fqmc::jagp_local_energy_adjoints_aabb(no, nv, &rt.t()(0,0), &tc(0,0), &rtc(0,0), &cur_im.t()(0,0), &cur_uum(0,0),
                                                        Jsi, Jsa, Jsj, Jsb, Jdia, Jdij, Jdib, Jdaj, Jdab, Jdjb,
                                                        &Tiajb.at(0), &ad_energy, &vrt(0,0), &vtc(0,0), &vt(0,0), &vartc(0,0),
                                                        &_Ksi.at(0), &_Ksa.at(0), &_Ksj.at(0), &_Ksb.at(0),
                                                        &_Kdia.at(0), &_Kdij.at(0), &_Kdib.at(0), &_Kdaj.at(0), &_Kdab.at(0), &_Kdjb.at(0));
        } else {
          formic::fqmc::agp_local_energy_adjoints_aabb(no, nv, &rt.t()(0,0), &tc(0,0), &rtc(0,0), &cur_im.t()(0,0), &cur_uum(0,0),
                                                       &Tiajb.at(0), &ad_energy, &vrt(0,0), &vtc(0,0), &vt(0,0), &vartc(0,0));
        }
        //formic::stop_timer("der_hc_autodiff_fortran");

        // add energy contribution
        local_e += ad_energy;

        // transpose the vtc matrix so it has dimensions (nv,no) (alpha,alpha)
        vtc.tip();

        // transpose the vartc matrix so it has dimensions (no,no) (beta,alpha)
        vartc.tip();

        // get derivatives resulting from vtc as the adjoint to rtheta
        this->rtheta_adjoint_to_der_vec(no, nv, ao2s, bo2s, au2s, bu2s, vtc, cur_im, cur_urm, der_out_a);

        // get derivatives resulting from vrt as the adjoint to theta*c
        this->rtheta_adjoint_to_der_vec(no, nv, bo2s, ao2s, bu2s, au2s, vrt.t(), cur_im.t(), cur_ucm.t(), der_out_b);

        // get derivatives resulting from vartc as the adjoint to theta (theta is the inverse pairing matrix)
        this->ipm_adjoint_to_der_vec(no, nv, ao2s, bo2s, vartc, cur_im.t(), der_out_a);

        // get derivatives resulting from vt as the adjoint to the unoccupied pairing matrix
        this->upm_adjoint_to_der_vec(no, nv, au2s, bu2s, vt, der_out_a);

        // change the sign of vt for the remaining two terms
        vt *= -formic::unity(S());

        // get derivatives resulting from vt*c as the adjoint to rtheta
        this->rtheta_adjoint_to_der_vec(no, nv, ao2s, bo2s, au2s, bu2s, vt * cur_ucm.t(), cur_im, cur_urm, der_out_a);

        // get derivatives resulting from vt* rtheta as the adjoint to the unoccupied column matrix
        this->urm_adjoint_to_der_vec(no, nv, bu2s, ao2s, vt.t() * rt, der_out_b);

        // get derivatives resulting from jastrows
        if (_apj) {
          //formic::start_timer("der_hc_autodiff_unpack_jf");
          S * aa_pos_a;
          S * bb_pos_a;
          S * ab_pos_a;
          _apj->get_aa_bb_ab_position_pointers(der_out_a + this->jf_pos(), aa_pos_a, bb_pos_a, ab_pos_a);
          S * aa_pos_b;
          S * bb_pos_b;
          S * ab_pos_b;
          _apj->get_aa_bb_ab_position_pointers(der_out_b + this->jf_pos(), aa_pos_b, bb_pos_b, ab_pos_b);
          _apj->intermed_adj_to_der_vec_ab(&ao2s[0], &bo2s[0], &au2s[0], &bu2s[0],
                                           Jsi, Jsa, Jsj, Jsb, Jdia, Jdij, Jdib, Jdaj, Jdab, Jdjb,
                                           &_Ksi.at(0), &_Ksa.at(0), &_Ksj.at(0), &_Ksb.at(0),
                                           &_Kdia.at(0), &_Kdij.at(0), &_Kdib.at(0), &_Kdaj.at(0), &_Kdab.at(0), &_Kdjb.at(0),
                                           aa_pos_a, ab_pos_a, aa_pos_b, ab_pos_b);
          //formic::stop_timer("der_hc_autodiff_unpack_jf");
        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Returns the Hamiltonian contraction with a derivative ratio vector.  Note that the first
      ///        element of the vector corresponds to the no-derivative term.
      ///        This version of the der_hc function uses the adjoint differentiation code.
      ///
      /// \param[in,out]  dfc      the derivative function coefficients.  size == 1+this->n_der()
      ///
      /// \return  sum_i <n|H|Psi_i> dfc_i / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S der_hc_autodiff(const S * const dfc) {

        // start a timer
        formic::start_timer("der_hc_autodiff");

        // make sure the object is initialized
        this->check_initialized("RGemPow::der_hc_autodiff");

        // get derivative ratios
        formic::Vector<S,int> drats(this->n_der(), formic::zero(S()));
        this->der_ratios(&drats.at(0));

        // get the local energy and its derivatives
        S local_e;
        formic::Vector<S,int> le_ders(this->n_der(), formic::zero(S()));
        this->local_energy_derivatives(local_e, &le_ders.at(0));

        // compute the contraction
        S retval = dfc[0] * local_e;
        for (int p = 0; p < le_ders.size(); p++)
          retval += ( local_e * drats[p] + le_ders[p] ) * dfc[ 1 + p ];
        //for (int p = this->pm_pos(); p < this->pm_pos()+_norb*_norb; p++)
        //  retval += ( local_e * drats[p] + le_ders[p] ) * dfc[ 1 + p ];
        //if (_apj)
        //  for (int p = this->jf_pos(); p < this->jf_pos()+_apj->n_der(); p++)
        //    retval += ( local_e * drats[p] + le_ders[p] ) * dfc[ 1 + p ];

        // stop the timer
        formic::stop_timer("der_hc_autodiff");

        // return the overall contraction
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the derivatives of the local energy w.r.t. the aa one-electron integral
      ///         and aaaa two-electron integral elements for excitation terms.
      ///         Non-excitation terms are handled separately.
      ///
      /// \param[in]      no       the number of occupied orbitals
      /// \param[in]      nv       the number of unoccupied orbitals
      /// \param[in]      ao2s     vector of occupied alpha sites
      /// \param[in]      bo2s     vector of occupied beta sites
      /// \param[in]      au2s     vector of unoccupied alpha sites
      /// \param[in]      bu2s     vector of unoccupied beta sites
      /// \param[in]      cur_im   the current inverse of the occupied pairing matrix
      /// \param[in]      cur_pm   the current all-site pairing matrix
      /// \param[in]      Jsi      occ jastrow intermediate
      /// \param[in]      Jsa      unocc jastrow intermediate
      /// \param[in]      Jdia     occ-unocc jastrow intermediate
      /// \param[in]      Jdai     unocc-occ jastrow intermediate (just the transpose of the Jdia matrix)
      /// \param[in]      Jdij     occ-occ jastrow intermediate
      /// \param[in]      Jdab     unocc-unocc jastrow intermediate
      /// \param[in]      weight   weight to apply to the derivatives
      /// \param[in,out]  mei_adj  cumulative adjoint derivatives for modiefied one electron integrals
      /// \param[in,out]  tei_adj  cumulative adjoint derivatives for two electron integrals
      /// \param[in,out]  local_e  the local energy
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void loc_en_integral_der_aa_aaaa(const int no,
                                       const int nv,
                                       const std::vector<int> & ao2s,
                                       const std::vector<int> & bo2s,
                                       const std::vector<int> & au2s,
                                       const std::vector<int> & bu2s,
                                       const ConstMat & cur_im,
                                       const ConstMat & cur_pm,
                                       const S * const  Jsi,
                                       const S * const  Jsa,
                                       const S * const Jdia,
                                       const S * const Jdai,
                                       const S * const Jdij,
                                       const S * const Jdab,
                                       const double weight,
                                             S * const mei_adj,
                                             S * const tei_adj,
                                       S & local_e) {

        // prepare bare aa and aaaa integrals
        this->prep_local_energy_ints_aa_aaaa(ao2s, bo2s, au2s, bu2s, _Tia_aa, _Tiajb_aaaa);

        // get the unoccupied row elements of the current pairing matrix
        //const ConstMat current_urm = cur_pm.slice(au2s, bo2s);

        // prepare contraction of new row matrix with inverse of occupied pairing matrix
        const ConstMat rt = cur_pm.slice(au2s, bo2s) * cur_im;
        //                   (new row mat)           * (inverse)
        //const ConstMat rt = current_urm * cur_im;

        // compute energy and adjoints
        S ad_energy = formic::zero(S());
        if (_apj) {
          formic::fqmc::jagp_local_energy_integral_der_aa_aaaa(no, nv, &rt.at(0,0), &rt.t().at(0,0),
                                                               Jsi, Jsa, Jdia, Jdai, Jdij, Jdab,
                                                               &_Tia_aa.at(0), &_Tiajb_aaaa.at(0), &ad_energy,
                                                               &_Tia_bb.at(0), &_Tiajb_bbbb.at(0));
        } else {
          throw formic::Exception("fomric::RGemPow::loc_en_integral_der_aa_aaaa currently only implemented for jagp");
        }

        // add energy contribution
        local_e += ad_energy;

        // add ov and ovov derivative adjoints to the integral derivatives
        this->add_local_energy_int_der_aa_aaaa(weight, ao2s, bo2s, au2s, bu2s, _Tia_bb, _Tiajb_bbbb, mei_adj, tei_adj);

      }

      void add_s_sqr_pentalty_to_aabb_ints(const int no,
                                           const int nv,
                                           const std::vector<int> & ao2s,
                                           const std::vector<int> & bo2s,
                                           const std::vector<int> & au2s,
                                           const std::vector<int> & bu2s,
                                                 S * const Vaabb) {

        if ( _s_sqr_penalty != 0.0 ) {
          for (int b = 0; b < nv; b++)
          for (int j = 0; j < no; j++)
          for (int a = 0; a < nv; a++)
          for (int i = 0; i < no; i++) {
            const double ib_equal = ( ao2s[i] - bu2s[b] == 0 ? 1.0 : 0.0 );
            const double ja_equal = ( bo2s[j] - au2s[a] == 0 ? 1.0 : 0.0 );
            const double ij_diff  = ( ao2s[i] - bo2s[j] == 0 ? 0.0 : 1.0 );
            Vaabb[ i + no * ( a + nv * ( j + no * b ) ) ] -= _s_sqr_penalty * ib_equal * ja_equal * ij_diff * formic::unity(S());
            //if ( _s_sqr_penalty * ib_equal * ja_equal * ij_diff != 0.0 )
            //  formic::of << boost::format("S^2 penalty added to %2i %2i %2i %2i") % ao2s[i] % au2s[a] % bo2s[j] % bu2s[b] << std::endl;
          }
        }
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the aabb two-electron integral contributions to the local energy and its
      ///         derivatives w.r.t. pairing matrix and jastrow variables.
      ///
      /// \param[in]      no         the number of occupied orbitals
      /// \param[in]      nv         the number of unoccupied orbitals
      /// \param[in]      ao2s       vector of occupied alpha sites
      /// \param[in]      bo2s       vector of occupied beta sites
      /// \param[in]      au2s       vector of unoccupied alpha sites
      /// \param[in]      bu2s       vector of unoccupied beta sites
      /// \param[in]      cur_im     the current inverse of the occupied pairing matrix
      /// \param[in]      cur_pm     the current all-site pairing matrix
      /// \param[in]      Jsi        jastrow single-index intermediate for index i
      /// \param[in]      Jsa        jastrow single-index intermediate for index a
      /// \param[in]      Jsj        jastrow single-index intermediate for index j
      /// \param[in]      Jsb        jastrow single-index intermediate for index b
      /// \param[in]      Jdia       jastrow double-index intermediate for indices i and a
      /// \param[in]      Jdij       jastrow double-index intermediate for indices i and j
      /// \param[in]      Jdib       jastrow double-index intermediate for indices i and b
      /// \param[in]      Jdaj       jastrow double-index intermediate for indices a and j
      /// \param[in]      Jdab       jastrow double-index intermediate for indices a and b
      /// \param[in]      Jdjb       jastrow double-index intermediate for indices j and b
      /// \param[in]      weight   weight to apply to the derivatives
      /// \param[in,out]  tei_adj  cumulative adjoint derivatives for two electron integrals
      /// \param[in,out]  local_e  the local energy
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void loc_en_integral_der_aabb(const int no,
                                    const int nv,
                                    const std::vector<int> & ao2s,
                                    const std::vector<int> & bo2s,
                                    const std::vector<int> & au2s,
                                    const std::vector<int> & bu2s,
                                    const ConstMat & cur_im,
                                    const ConstMat & cur_pm,
                                    const S * const  Jsi,
                                    const S * const  Jsa,
                                    const S * const  Jsj,
                                    const S * const  Jsb,
                                    const S * const Jdia,
                                    const S * const Jdij,
                                    const S * const Jdib,
                                    const S * const Jdaj,
                                    const S * const Jdab,
                                    const S * const Jdjb,
                                    const double weight,
                                          S * const tei_adj,
                                    S & local_e) {

        // prepare bare aabb integrals
        formic::tei_get_ovov_aabb_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), _ham.dense_tei_ptr(), &_Tiajb_aaaa.at(0));

        // add S^2 penalty if necessary
        this->add_s_sqr_pentalty_to_aabb_ints(no, nv, ao2s, bo2s, au2s, bu2s, &_Tiajb_aaaa.at(0));

        // get the unoccupied row elements of the current pairing matrix
        const ConstMat cur_urm = cur_pm.slice(au2s, bo2s);

        // get the unoccupied column elements of the current pairing matrix
        const ConstMat cur_ucm = cur_pm.slice(ao2s, bu2s);

        // get the unocc-unocc pairing matrix
        const ConstMat cur_uum = cur_pm.slice(au2s, bu2s);

        // prepare rt(a,i) = contraction of new row matrix with inverse of occupied pairing matrix
        const ConstMat rt = cur_urm * cur_im;

        // prepare tc(j,b) = contraction of new column matrix with inverse of occupied pairing matrix
        const ConstMat tc = cur_im * cur_ucm;

        // prepare rtc(a,b) = contraction of new row, inverse, and new column matrices
        const ConstMat rtc = cur_urm * tc;

        // compute energy and adjoints
        S ad_energy = formic::zero(S());
        if (_apj) {
          formic::fqmc::jagp_local_energy_integral_der_aabb(no, nv, &rt.t().at(0,0), &tc.at(0,0), &rtc.at(0,0), &cur_im.t().at(0,0), &cur_uum.at(0,0),
                                                            Jsi, Jsa, Jsj, Jsb, Jdia, Jdij, Jdib, Jdaj, Jdab, Jdjb,
                                                            &_Tiajb_aaaa.at(0), &ad_energy, &_Tiajb_bbbb.at(0));
        } else {
          throw formic::Exception("fomric::RGemPow::loc_en_integral_der_aabb currently only implemented for jagp");
        }

        // add energy contribution
        local_e += ad_energy;

        // add integral adjoint derivative contributions to total
        const int nt = no + nv;
        for (int b = 0; b < nv; b++)
        for (int j = 0; j < no; j++)
        for (int a = 0; a < nv; a++)
        for (int i = 0; i < no; i++) {
          S & integral = tei_adj[ ao2s[i] + nt * ( au2s[a] + nt * ( bo2s[j] + nt * bu2s[b] ) ) ];
          integral += weight * _Tiajb_bbbb[ i + no * ( a + nv * ( j + no * b ) ) ];
        }
        //const int nt = no + nv;
        //for (int b = 0; b < nv; b++)
        //for (int j = 0; j < no; j++)
        //for (int a = 0; a < nv; a++)
        //for (int i = 0; i < no; i++) {
        //  const S integral = tei_dense[ ao2s[i] + nt * ( au2s[a] + nt * ( bo2s[j] + nt * bu2s[b] ) ) ];
        //  tei_slice[ i + no * ( a + nv * ( j + no * b ) ) ] = integral;
        //}
        //formic::tei_get_ovov_aabb_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), _ham.dense_tei_ptr(), &_Tiajb_aaaa.at(0));

      }

    // protected member functions
    protected:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief updates the coefficient object's variables
      ///
      /// \param[in]      dfc      the derivative function coefficients.  size == this->n_der()
      /// \param[in]      alpha    the constant to scale the update by
      ///
      ///  |Psi(x_1, x_2, ...)> + sum_{i=1} dfc_i |dPsi/dx_i>
      ///     --->  |Psi(x_1 + alpha * dfc_1, x_2 + alpha * dfc_2, ...)>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_variables_detail(const S * const dfc, const S alpha) {

        // add pairing matrix derivatives
        if (!_fixed_pm) {
          formic::xaxpy(_norb*_norb, alpha, dfc + this->pm_pos(), 1, _pm, 1);
          //const S * const pm_der = dfc + this->pm_pos();
          //for (int i = 0; i < _norb*_norb; i++)
          //  _pm[i] += alpha * pm_der[i];
        }

        //// symmetrize pairing matrix
        //for (int i = 0; i < _norb; i++)
        //  for (int j = 0; j < i; j++)
        //    _pm[_norb*i+j] = _pm[_norb*j+i];

        //// scale pairing matrix so that the maximum-magnitude element is one
        //if (!_fixed_pm) {
        //  S max_e = formic::zero(S());
        //  for (int i = 0; i < _norb*_norb; i++)
        //    if ( std::abs(_pm[i]) > std::abs(max_e) )
        //      max_e = _pm[i];
        //  formic::xscal(_norb*_norb, 1.0 / max_e, _pm, 1);
        //}

        // update jastrow variables
        if (_apj)
          _apj->update_variables(dfc + this->jf_pos(), alpha);

        // clear the exact sampling vectors
        this->clear_exact_sample_vecs();

      }

    // public member functions
    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief constructs the coefficient object by reading names from an archive and finding the
      ///        pairing matrix in the supplied wavefunction
      ///
      /// \param[in]      userinp  user input options
      /// \param[in]      ham      the Hamiltonian
      /// \param[in,out]  wfn      the wavefunction holding the pairing matrix data
      /// \param[in,out]  arch     an archive from which to read names
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      RGemPow(const formic::InputBase & userinp,
              formic::Hamiltonian<S> & ham,
              formic::Wavefunction<S> & wfn,
              formic::Archive & arch,
              const bool have_jastrow)
        : _norb(userinp.get<int>("nao")),
          _nap(userinp.get<int>("nalpha")),
          _opm(userinp.get<int>("nalpha") * userinp.get<int>("nalpha"), formic::zero(S())),
          _oim(userinp.get<int>("nalpha") * userinp.get<int>("nalpha"), formic::zero(S())),
          _oic(userinp.get<int>("nalpha") * userinp.get<int>("nalpha"), formic::zero(S())),
          _swork( std::max(5, userinp.get<int>("nalpha")) * userinp.get<int>("nalpha"), formic::zero(S())),
          _iwork(3 * userinp.get<int>("nalpha"), 0),
          _ham(ham),
          _initialized(false),
          _s_sqr_penalty(userinp.get<double>("s_sqr_penalty")),
          _config(userinp),
          _apj(0),
          _der_hc_by_adj(userinp.get<bool>("der_hc_by_adj"))
      {

        // ensure there are equal alpha and beta electron numbers
        if ( userinp.get<int>("nalpha") != userinp.get<int>("nbeta") )
          throw formic::Exception("RGemPow coeff object assumes equal numbers of alpha and beta electrons");

        // read in the name to give the coefficient object
        formic::mpi::read_and_bcast(arch, this->_name, "Failed to read name in RGemPow constructor.");

        // read in name of pairing matrix
        std::string pmn;
        formic::mpi::read_and_bcast(arch, pmn, "Failed to read name of pairing matrix data object in RGemPow constructor.");

        // if requested, get the jastrow factors
        if (have_jastrow) {
          std::string jname;
          formic::mpi::read_and_bcast(arch, jname, "Failed to read name of jastrow factor coefficient object in RGemPow constructor.");
          formic::CoeffObj<S> * const jast_co_ptr = wfn.get_co(jname).get();
          if ( jast_co_ptr == 0 )
            throw formic::Exception("pointer for all-pair jastrow factor coefficient object \"%s\" was empty in constructor for RGemPow coeff object \"%s\"") % jname % this->_name;
          _apj = dynamic_cast<formic::fqmc::AllPairJast<S> *>(jast_co_ptr);
          if ( _apj == 0 )
            throw formic::Exception("casted pointer for all-pair jastrow factor coefficient object \"%s\" is invalid in constructor for RGemPow coeff object \"%s\"") % jname % this->_name;
        }

        // get a reference to the vector holding the pairing matrix
        std::vector<S> & pmv = wfn.get_wf(pmn)->data_vec();

        // ensure the pairing matrix is the correct size
        if ( pmv.size() !=  _norb * _norb )
          throw formic::Exception("pairing matrix \"%s\" is the wrong size in constructor for RGemPow coeff object \"%s\"") % pmn % this->_name;

        // set the pairing matrix pointer
        assert( pmv.size() > 0 );
        _pm = &pmv.at(0);

        // get whether the pairing matrix is fixed
        _fixed_pm = wfn.get_wf(pmn)->fixed();

        //// symmetrize pairing matrix
        //for (int i = 0; i < _norb; i++)
        //  for (int j = 0; j < i; j++)
        //    _pm[_norb*i+j] = _pm[_norb*j+i];

        // initialize tensors for local energy and derivative contraction
        {

          const int no = _nap;
          const int nv = _norb - _nap;

          _Tia_aa.assign(             no*nv, formic::zero(S()));
          _Tia_bb.assign(             no*nv, formic::zero(S()));
          _Tiajb_aaaa.assign(   no*nv*no*nv, formic::zero(S()));
          _Tiajb_bbbb.assign(   no*nv*no*nv, formic::zero(S()));
          _Tiajb_aabb.assign(   no*nv*no*nv, formic::zero(S()));
          _Bkl.assign(               no*no, formic::zero(S()));
          _Blk.assign(               no*no, formic::zero(S()));
          _Via.assign(               no*nv, formic::zero(S()));
          _Ujb.assign(               no*nv, formic::zero(S()));
          _Pab.assign(               nv*nv, formic::zero(S()));
          _VBia.assign(               no*nv, formic::zero(S()));
          _BUjb.assign(               no*nv, formic::zero(S()));
          _TVBjb_aaaa.assign(                no*nv, formic::zero(S()));
          _TVBjb_bbbb.assign(                no*nv, formic::zero(S()));
          _TVBjb_aabb.assign(                no*nv, formic::zero(S()));
          _TVBaj.assign(                no*nv, formic::zero(S()));
          _TUBaj.assign(                no*nv, formic::zero(S()));
          _VBUab.assign(              nv*nv, formic::zero(S()));
          _TBab.assign(               nv*nv, formic::zero(S()));

          _Ckl.assign(               no*no, formic::zero(S()));
          _Clk.assign(               no*no, formic::zero(S()));
          _Zia.assign(               no*nv, formic::zero(S()));
          _Yjb.assign(               no*nv, formic::zero(S()));
          _Qab.assign(               nv*nv, formic::zero(S()));
          _ZBia.assign(               no*nv, formic::zero(S()));
          _VCia.assign(               no*nv, formic::zero(S()));
          _BYjb.assign(               no*nv, formic::zero(S()));
          _CUjb.assign(               no*nv, formic::zero(S()));
          _dVBUab.assign(               nv*nv, formic::zero(S()));
          _TCab.assign(               nv*nv, formic::zero(S()));
          _TBUia.assign(               no*nv, formic::zero(S()));

          _RTA.assign(  no*nv, formic::zero(S()));
          _Ksi.assign(     no, formic::zero(S()));
          _Ksj.assign(     no, formic::zero(S()));
          _Ksa.assign(     nv, formic::zero(S()));
          _Ksb.assign(     nv, formic::zero(S()));
          _Kdia.assign( no*nv, formic::zero(S()));
          _Kdij.assign( no*no, formic::zero(S()));
          _Kdib.assign( no*nv, formic::zero(S()));
          _Kdaj.assign( no*nv, formic::zero(S()));
          _Kdab.assign( nv*nv, formic::zero(S()));
          _Kdjb.assign( no*nv, formic::zero(S()));

        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the name of the child class
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      std::string class_name() const { return "RGemPow"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of derivatives, which is the number of variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      int n_der() const {
        int retval = ( _fixed_pm ? 0 : _norb * _norb );
        if (_apj)
          retval += _apj->n_der();
        return retval;
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of variables (counting the "fixed" ones)
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      int n_var() const {
        int retval = _norb * _norb;
        if (_apj)
          retval += _apj->n_var();
        return retval;
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Initializes the object for the given lattice configuration.
      ///
      /// \param[in]      config  the lattice configuration we want the object initialized for
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void initialize(const formic::Configuration & config) {

        // save the configuration
        _config = config;

        // construct the occupied pairing matrix
        formic::fqmc::rgp_build_occupied_matrix(_nap, _norb, _pm,
                                                &_config.a_occ_to_site().at(0),
                                                &_config.b_occ_to_site().at(0),
                                                &_opm.at(0));

        // compute the inverse of the occupied pairing matrix
        formic::xcopy(_opm.size(), &_opm.at(0), 1, &_oim.at(0), 1);
        S det_value = formic::unity(S());
        formic::matrix_inverse_lu(_nap, det_value, &_oim.at(0), &_swork.at(0), &_iwork.at(0));

        // initialize the all pair jastrows
        if (_apj)
          _apj->initialize(config);

        // set the current coefficient
        _curr_coeff = this->get_coeff(config);

        // set the initialized flag to true
        _initialized = true;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Updates the object so that the current configuration becomes the new configuration,
      ///        i.e. current = current + diff
      ///
      /// \param[in]      diff    an object describing the change to the configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update(const formic::ConfigDiff & diff) {

        // check that the object is initialized
        this->check_initialized("RGemPow::update");

        // check that only one or two electrons moved
        if ( diff.newly_full().size() != 1 && diff.newly_full().size() != 2 )
          throw formic::Exception("there should be 1 or 2 newly full sites in RGemPow::update");
        if ( diff.newly_empty().size() != 1 && diff.newly_empty().size() != 2 )
          throw formic::Exception("there should be 1 or 2 newly empty sites in RGemPow::update");

        // check that empty and full are the same length
        if ( diff.newly_empty().size() != diff.newly_full().size() )
          throw formic::Exception("there should be equal numbers of newly full and empty sites in RGemPow::update");

        // check that the old/new occupied sites are of the same spin
        for (int i = 0; i < diff.newly_full().size(); i++)
          if ( diff.newly_full().at(i) % 2 != diff.newly_empty().at(i) % 2 )
            throw formic::Exception("new/old site are of different spins in RGemPow::update");

        // check that the configuration change is possible
        formic::check_config_update(diff, _config, "RGemPow::update");
        //this->check_config_diff(diff, "RGemPow::update");

        // update the current coefficient
        _curr_coeff *= this->get_ratio(diff);

        // update the all pair jastrows
        if (_apj)
          _apj->update(diff);

        bool config_updated = false;

        // alpha move
        if ( diff.newly_full().size() == 1 && diff.newly_full().at(0) % 2 == 0 ) {

          // update the occupied pairing matrix and its inverse
          formic::fqmc::rgp_row_update(_nap,
                                       _norb,
                                       diff.newly_full().at(0)/2,
                                       diff.newly_empty().at(0)/2,
                                       _pm,
                                       &_config.a_site_to_occ().at(0),
                                       &_config.a_occ_to_site().at(0),
                                       &_config.b_occ_to_site().at(0),
                                       &_opm.at(0),
                                       &_oim.at(0),
                                       &_swork.at(0*_nap),
                                       &_swork.at(1*_nap),
                                       &_swork.at(2*_nap));

        // beta move
        } else if ( diff.newly_full().size() == 1 && diff.newly_full().at(0) % 2 == 1 ) {

          // update the occupied pairing matrix and its inverse
          formic::fqmc::rgp_col_update(_nap,
                                       _norb,
                                       diff.newly_full().at(0)/2,
                                       diff.newly_empty().at(0)/2,
                                       _pm,
                                       &_config.b_site_to_occ().at(0),
                                       &_config.b_occ_to_site().at(0),
                                       &_config.a_occ_to_site().at(0),
                                       &_opm.at(0),
                                       &_oim.at(0),
                                       &_swork.at(0*_nap),
                                       &_swork.at(1*_nap),
                                       &_swork.at(2*_nap));

        // double moves
        } else {

          _config.update(diff);
          this->initialize(_config);
          config_updated = true;

        }

        // finally, update the configuration
        if (!config_updated)
          _config.update(diff);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      formic::MantExp<S> get_coeff() {

        // check that the object is initialized
        this->check_initialized("RGemPow::get_coeff()");

        return _curr_coeff;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for an arbitrary configuration
      ///
      /// \param[in]      config  the lattice configuration for which we want the coefficient
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      formic::MantExp<S> get_coeff(const formic::Configuration & config) {

        // build the occupied pairing matrix for the supplied configuration
        formic::fqmc::rgp_build_occupied_matrix(_nap, _norb, _pm, &config.a_occ_to_site().at(0), &config.b_occ_to_site().at(0), &_swork.at(0));

        // compute the determinant of the pairing matrix
        formic::MantExp<S> retval = formic::matrix_determinant_lu(_nap, &_swork.at(0), &_iwork.at(0));

        // add jastrow factor contribution
        if (_apj)
          retval = retval * _apj->get_coeff(config);

        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the ratio of the coefficient at a new configuration (current + diff) to
      ///        that at the current configuration
      ///
      /// \param[in]      diff     the differece between the new and current configurations
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S get_ratio(const formic::ConfigDiff & diff) {

        // check that the object is initialized
        this->check_initialized("RGemPow::get_ratio");

        // check that empty and full lists are the same size
        if ( diff.newly_full().size() != diff.newly_empty().size() )
          throw formic::Exception("different sized newly empty and newly full lists in RGemPow::get_ratio");

        // check that the configuration change is possible
        formic::check_config_update(diff, _config, "RGemPow::get_ratio");
        //this->check_config_diff(diff, "RGemPow::get_ratio");

        // get the geminal power ratio
        S retval = formic::fqmc::rgp_det_ratio(_nap,
                                               _norb,
                                               2*diff.newly_empty().size(),
                                               ( diff.newly_full().size()  > 0 ? &diff.newly_full().at(0)  : 0 ),
                                               ( diff.newly_empty().size() > 0 ? &diff.newly_empty().at(0) : 0 ),
                                               _pm,
                                               &_oim.at(0),
                                               &_config.a_occ_to_site().at(0),
                                               &_config.b_occ_to_site().at(0),
                                               &_config.a_site_to_occ().at(0),
                                               &_config.b_site_to_occ().at(0),
                                               &_swork.at(0*_nap),
                                               &_swork.at(1*_nap),
                                               &_swork.at(2*_nap));

        // get the jastrow ratios
        if (_apj)
          retval = retval * _apj->get_ratio(diff);

        // return the overall ratio
        return retval;

      }

//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns part of the local energy for the current configuration
//      ///        (for development and debugging)
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      S partial_local_energy() {
//
//        // make sure the object is initialized
//        this->check_initialized("RGemPow::local_energy");
//
//        // initialize local energy
//        S retval = formic::zero(S());
//
//        // get number of ocupied and unoccuied sites
//        const int no = _nap;
//        const int nv = _norb - _nap;
//
//        // get references to occupation arrays
//        const std::vector<int> & ao2s = _config.a_occ_to_site();
//        const std::vector<int> & bo2s = _config.b_occ_to_site();
//        const std::vector<int> & au2s = _config.a_unocc_to_site();
//        const std::vector<int> & bu2s = _config.b_unocc_to_site();
//
//        // prepare inverse matrix and transpose
//        for (int j = 0; j < no; j++)
//        for (int i = 0; i < no; i++) {
//          _Bkl[i+no*j] = _oim[i*no+j];
//          _Blk[j+no*i] = _oim[i*no+j];
//        }
//
//        // prepare new row, new column, and unoccupied pairing matrices
//        this->prep_row_col_unocc_mat(_pm, _Via, _Ujb, _Pab);
//
//        // prepare bare aa and aaaa integrals
//        this->prep_local_energy_ints_aa_aaaa(ao2s, bo2s, au2s, bu2s, _Tia_aa, _Tiajb_aaaa);
//        std::fill(_Tia_aa.begin(), _Tia_aa.end(), formic::zero(S()));
//        std::fill(_Tiajb_aaaa.begin(), _Tiajb_aaaa.end(), formic::zero(S()));
//
//        // transform the aa and aaaa integrals by the Jastrow factor ratios
//        if (_apj)
//          _apj->transform_ints_aa_aaaa(&_Tia_aa.at(0), &_Tiajb_aaaa.at(0));
//
//        // compute aa and aaaa one and two electron excitation contribution
//        // (note that this also prepares VBia, which will be used for aabb below)
//        formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_aa[0], &_Tiajb_aaaa[0], &_Bkl[0], &_Via[0],
//                                      &_VBia[0], &_TVBjb_aaaa[0], &_TVBaj[0], &retval);
//
//        // prepare bare bb and bbbb integrals
//        this->prep_local_energy_ints_aa_aaaa(bo2s, ao2s, bu2s, au2s, _Tia_bb, _Tiajb_bbbb);
//        std::fill(_Tia_bb.begin(), _Tia_bb.end(), formic::zero(S()));
//        std::fill(_Tiajb_bbbb.begin(), _Tiajb_bbbb.end(), formic::zero(S()));
//
//        // transform the bb and bbbb integrals by the Jastrow factor ratios
//        if (_apj)
//          _apj->transform_ints_bb_bbbb(&_Tia_bb.at(0), &_Tiajb_bbbb.at(0));
//
//        // compute bb and bbbb one and two electron excitation contribution
//        // (note that this also prepares BUjb, which will be used for aabb below)
//        formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_bb[0], &_Tiajb_bbbb[0], &_Blk[0], &_Ujb[0],
//                                      &_BUjb[0], &_TVBjb_bbbb[0], &_TUBaj[0], &retval);
//
////        // prepare bare aabb integrals
////        formic::tei_get_ovov_aabb_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), _ham.dense_tei_ptr(), &_Tiajb_aabb.at(0));
////
////        // transform the aabb integrals by the Jastrow factor ratios
////        if (_apj)
////          _apj->transform_ints_aabb(&_Tiajb_aabb.at(0));
////
////        // compute aabb contribution
////        formic::fqmc::rgp_aabb_energy(no, nv, &_Tiajb_aabb[0], &_Via[0], &_Ujb[0], &_Pab[0],
////                                      &_Bkl[0], &_Blk[0], &_VBia[0], &_BUjb[0],
////                                      &_VBUab[0], &_TBab[0], &_TVBjb_aabb[0], &retval);
//
//        // add non-excitation terms
//        retval += this->non_excite_energy_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s);
//        retval += this->non_excite_energy_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s);
//
//        // add nuclear repulsion term
//        retval += _ham.e_cor_nuc() * formic::unity(S());
//
//        // return the local energy
//        return retval;
//
//      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the local energy for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S local_energy() {

        // make sure the object is initialized
        this->check_initialized("RGemPow::local_energy");

        // initialize local energy
        S retval = formic::zero(S());

        // get number of ocupied and unoccuied sites
        const int no = _nap;
        const int nv = _norb - _nap;

        // get references to occupation arrays
        const std::vector<int> & ao2s = _config.a_occ_to_site();
        const std::vector<int> & bo2s = _config.b_occ_to_site();
        const std::vector<int> & au2s = _config.a_unocc_to_site();
        const std::vector<int> & bu2s = _config.b_unocc_to_site();

        // prepare inverse matrix and transpose
        for (int j = 0; j < no; j++)
        for (int i = 0; i < no; i++) {
          _Bkl[i+no*j] = _oim[i*no+j];
          _Blk[j+no*i] = _oim[i*no+j];
        }

        // prepare new row, new column, and unoccupied pairing matrices
        this->prep_row_col_unocc_mat(_pm, _Via, _Ujb, _Pab);

        // prepare bare aa and aaaa integrals
        this->prep_local_energy_ints_aa_aaaa(ao2s, bo2s, au2s, bu2s, _Tia_aa, _Tiajb_aaaa);

        // transform the aa and aaaa integrals by the Jastrow factor ratios
        if (_apj)
          _apj->transform_ints_aa_aaaa(&_Tia_aa.at(0), &_Tiajb_aaaa.at(0));

        // compute aa and aaaa one and two electron excitation contribution
        // (note that this also prepares VBia, which will be used for aabb below)
        formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_aa[0], &_Tiajb_aaaa[0], &_Bkl[0], &_Via[0],
                                      &_VBia[0], &_TVBjb_aaaa[0], &_TVBaj[0], &retval);

        // prepare bare bb and bbbb integrals
        this->prep_local_energy_ints_aa_aaaa(bo2s, ao2s, bu2s, au2s, _Tia_bb, _Tiajb_bbbb);

        // transform the bb and bbbb integrals by the Jastrow factor ratios
        if (_apj)
          _apj->transform_ints_bb_bbbb(&_Tia_bb.at(0), &_Tiajb_bbbb.at(0));

        // compute bb and bbbb one and two electron excitation contribution
        // (note that this also prepares BUjb, which will be used for aabb below)
        formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_bb[0], &_Tiajb_bbbb[0], &_Blk[0], &_Ujb[0],
                                      &_BUjb[0], &_TVBjb_bbbb[0], &_TUBaj[0], &retval);

        // prepare bare aabb integrals
        formic::tei_get_ovov_aabb_slice(no, nv, &ao2s.at(0), &bo2s.at(0), &au2s.at(0), &bu2s.at(0), _ham.dense_tei_ptr(), &_Tiajb_aabb.at(0));

        // add S^2 penalty if necessary
        this->add_s_sqr_pentalty_to_aabb_ints(no, nv, ao2s, bo2s, au2s, bu2s, &_Tiajb_aabb.at(0));

        // transform the aabb integrals by the Jastrow factor ratios
        if (_apj)
          _apj->transform_ints_aabb(&_Tiajb_aabb.at(0));

        // compute aabb contribution
        formic::fqmc::rgp_aabb_energy(no, nv, &_Tiajb_aabb[0], &_Via[0], &_Ujb[0], &_Pab[0],
                                      &_Bkl[0], &_Blk[0], &_VBia[0], &_BUjb[0],
                                      &_VBUab[0], &_TBab[0], &_TVBjb_aabb[0], &retval);

        // add non-excitation terms
        retval += this->non_excite_energy_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s);
        retval += this->non_excite_energy_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s);

        // add nuclear repulsion term
        retval += _ham.e_cor_nuc() * formic::unity(S());

        // return the local energy
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the diagonal contribution to the local energy for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S local_diag_energy() {

        // initialize return value
        S retval = formic::zero(S());

        // one electron energy
        for (int i = 0; i < _nap; i++) {
          const int ai = _config.a_occ_to_site().at(i);
          const int bi = _config.b_occ_to_site().at(i);
          retval += _ham.get_mei(ai, ai);
          retval += _ham.get_mei(bi, bi);
        }

        // two electron coulomb energy
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _nap; j++) {
          const int ai = _config.a_occ_to_site().at(i);
          const int aj = _config.a_occ_to_site().at(j);
          const int bi = _config.b_occ_to_site().at(i);
          const int bj = _config.b_occ_to_site().at(j);
          retval += 0.5 * _ham.get_tei(ai, ai, aj, aj);
          retval += 0.5 * _ham.get_tei(bi, bi, bj, bj);
          retval += 1.0 * _ham.get_tei(ai, ai, bj, bj);
        }

        // two electron exchange energy
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _norb - _nap; j++) {
          const int ai = _config.a_occ_to_site().at(i);
          const int aj = _config.a_unocc_to_site().at(j);
          const int bi = _config.b_occ_to_site().at(i);
          const int bj = _config.b_unocc_to_site().at(j);
          retval += 0.5 * _ham.get_tei(ai, aj, aj, ai);
          retval += 0.5 * _ham.get_tei(bi, bj, bj, bi);
        }

        // add nuclear repulsion term
        retval += _ham.e_cor_nuc() * formic::unity(S());

        // return the diagonal contribution to the local energy
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief gives the derivative ratios for the current configuration
      ///
      /// \param[out]     ratios     on exit, the derivative ratios.
      ///                            for ign_fixed == true, size == this->n_der()
      ///                            for ign_fixed == false, size == this->n_var()
      /// \param[in]      ign_fixed  whether to ignore variables that are fixed
      ///
      ///   The ratio for variable i is     <n|Psi_i> / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void der_ratios(S * const ratios, const bool ign_fixed = true) {

        this->der_ratios_detail(true, ratios, ign_fixed);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Computes the local energy and its derivatives with respect to the wavefunction variables
      ///
      /// \param[out]     local_e  On exit, the local energy.
      /// \param[out]     le_ders  On exit, the derivatives of the local energy.  size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void local_energy_derivatives(S & local_e, S * const le_ders) {

        // make sure the object is initialized
        this->check_initialized("RGemPow::local_energy_derivatives");

        // get references to occupation arrays
        const std::vector<int> & ao2s = _config.a_occ_to_site();
        const std::vector<int> & bo2s = _config.b_occ_to_site();
        const std::vector<int> & au2s = _config.a_unocc_to_site();
        const std::vector<int> & bu2s = _config.b_unocc_to_site();

        // get space to hold the beta-first component of the derivatives
        formic::Vector<S,int> le_der_storage_b(this->n_der());

        // get pointers to the alpha-first and beta-first components
        S * const le_der_vec_a = le_ders;
        S * const le_der_vec_b = &le_der_storage_b.at(0);

        // initialize the derivatives to zero
        std::fill(le_der_vec_a, le_der_vec_a + this->n_der(), formic::zero(S()));
        std::fill(le_der_vec_b, le_der_vec_b + this->n_der(), formic::zero(S()));

        // get the current pairing matrix (transposing to switch from c to fortran indexing)
        const ConstMat current_pm = ConstMat(_norb, _norb, _pm).t();

        // get the inverse of the occupied pairing matrix (transposing to switch from c to fortran indexing)
        const ConstMat current_im = ConstMat(_nap, _nap, &_oim.at(0)).t();

        // initialize some local energy components
        S local_e_aaaa = formic::zero(S());
        S local_e_bbbb = formic::zero(S());
        S local_e_aabb = formic::zero(S());

        // compute aa/aaaa contribution to energy and derivatives
        if (_apj) {
          this->loc_en_der_aa_aaaa(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                   current_im, current_pm,
                                   _apj->Jsi_a_ptr(), _apj->Jsa_a_ptr(),
                                   _apj->Jdia_aa_ptr(), _apj->Jdai_aa_ptr(),
                                   _apj->Jdij_aa_ptr(), _apj->Jdab_aa_ptr(),
                                   _Tia_aa, _Tiajb_aaaa, local_e_aaaa, le_der_vec_a);
        } else {
          this->loc_en_der_aa_aaaa(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                   current_im, current_pm,
                                   NULL, NULL,
                                   NULL, NULL,
                                   NULL, NULL,
                                   _Tia_aa, _Tiajb_aaaa, local_e_aaaa, le_der_vec_a);
        }

        // compute bb/bbbb contribution to energy and derivatives
        // note that this is done by exchanging alpha and beta and using the alpha code
        if (_apj) {
          this->loc_en_der_aa_aaaa(_nap, _norb-_nap, bo2s, ao2s, bu2s, au2s,
                                   current_im.t(), current_pm.t(),
                                   _apj->Jsi_b_ptr(), _apj->Jsa_b_ptr(),
                                   _apj->Jdia_bb_ptr(), _apj->Jdai_bb_ptr(),
                                   _apj->Jdij_bb_ptr(), _apj->Jdab_bb_ptr(),
                                   _Tia_bb, _Tiajb_bbbb, local_e_bbbb, le_der_vec_b);
        } else {
          this->loc_en_der_aa_aaaa(_nap, _norb-_nap, bo2s, ao2s, bu2s, au2s,
                                   current_im.t(), current_pm.t(),
                                   NULL, NULL,
                                   NULL, NULL,
                                   NULL, NULL,
                                   _Tia_bb, _Tiajb_bbbb, local_e_bbbb, le_der_vec_b);
        }

        // compute aabb contribution to energy and derivatives
        if (_apj) {
          this->loc_en_der_aabb(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                current_im, current_pm,
                                _apj->Jsi_a_ptr(),
                                _apj->Jsa_a_ptr(),
                                _apj->Jsi_b_ptr(),
                                _apj->Jsa_b_ptr(),
                                _apj->Jdia_aa_ptr(),
                                _apj->Jdij_ab_ptr(),
                                _apj->Jdia_ab_ptr(),
                                _apj->Jdai_ab_ptr(),
                                _apj->Jdab_ab_ptr(),
                                _apj->Jdia_bb_ptr(),
                                _Tiajb_aabb, local_e_aabb, le_der_vec_a, le_der_vec_b);
        } else {
          this->loc_en_der_aabb(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                current_im, current_pm,
                                NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                                _Tiajb_aabb, local_e_aabb, le_der_vec_a, le_der_vec_b);
        }

        // combine beta-first and alpha-first components of the derivatives
        for (int i = 0; i < _norb; i++) // alpha
        for (int j = 0; j < _norb; j++) // beta
          le_der_vec_a[ this->pm_pos() + i * _norb + j ] += le_der_vec_b[ this->pm_pos() + j * _norb + i ];
        if (_apj)
          _apj->add_beta_derivs_to_alpha_derivs(&le_der_vec_a[this->jf_pos()], &le_der_vec_b[this->jf_pos()]);

        // combine the components of the local energy
        local_e =   local_e_aaaa
                  + local_e_bbbb
                  + local_e_aabb
                  + this->non_excite_energy_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s)
                  + this->non_excite_energy_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s)
                  + _ham.e_cor_nuc() * formic::unity(S());

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the orbital coeff matrix derivative ratios for the current configuration.
      ///         
      /// \param[in]      Cip      The orbital coefficient matrix.
      ///                          size == _norb*_norb
      /// \param[out]     dPsidCip On exit, the derivative ratios w.r.t. the orbital coefficient matrix.
      ///                          size == _norb*_norb
      ///
      ///         The orbital coefficient derivative ratio w.r.t. the (i,p) element of the orbital
      ///         coefficient matrix C is defined as  <dPsi/dC(i,p)|n> / <Psi|n>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void orb_coeff_der_ratios(const S * const Cip, S * const dPsidCip) {

        // make sure the object is initialized
        this->check_initialized("RGemPow::orb_coeff_der_ratios");

        // get the pairing matrix derivative ratios (here we skip the jastrows, defining them as fixed regardless of orbital rotation)
        Mat pm_dr(_norb, _norb, formic::zero(S()));
        this->der_ratios_detail(false, &pm_dr.at(0,0));

        // put the pairing matrix elements into a matrix (trasposing to switch from c to fortran indexing)
        const ConstMat pm_mat = ConstMat(_norb, _norb, _pm).t();

        // organize the orbital rotation elements into a matrix
        const ConstMat cip_mat(_norb, _norb, Cip);

        // use a matrix to initialize and organize the orbital rotation derivative ratios
        Mat ordr_mat(_norb, _norb, dPsidCip, formic::zero(S()));

        // compute the orbital rotation matrix derivative ratios
        ordr_mat += pm_dr * cip_mat * pm_mat.t();
        ordr_mat += pm_dr.t() * cip_mat * pm_mat;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Gets the local configuration's weighted contribution to the energy derivatives w.r.t.
      ///         the one and two electron integral elements.
      ///         
      /// \param[in]      weight   Weight to give the local contributions.
      /// \param[out]     local_e  On exit, the local energy.
      /// \param[in,out]  mei_adj  Cumulative derivative adjoints for modified one electron integrals.
      ///                          size == _norb*_norb
      /// \param[in,out]  tei_adj  Cumulative derivative adjoints for two electron integrals.
      ///                          size == _norb*_norb*_norb*_norb
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void local_energy_integral_der(const double weight, S & local_e, S * const mei_adj, S * const tei_adj) {

        // make sure the object is initialized
        this->check_initialized("RGemPow::local_energy_integral_der");

        // get references to occupation arrays
        const std::vector<int> & ao2s = _config.a_occ_to_site();
        const std::vector<int> & bo2s = _config.b_occ_to_site();
        const std::vector<int> & au2s = _config.a_unocc_to_site();
        const std::vector<int> & bu2s = _config.b_unocc_to_site();

        // get the current pairing matrix (transposing to switch from c to fortran indexing)
        const ConstMat current_pm = ConstMat(_norb, _norb, _pm).t();

        // get the inverse of the occupied pairing matrix (transposing to switch from c to fortran indexing)
        const ConstMat current_im = ConstMat(_nap, _nap, &_oim.at(0)).t();

        // initialize some local energy components
        S local_e_aaaa = formic::zero(S());
        S local_e_bbbb = formic::zero(S());
        S local_e_aabb = formic::zero(S());

        // compute aa/aaaa contribution to energy and derivatives
        if (_apj) {
          this->loc_en_integral_der_aa_aaaa(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                            current_im, current_pm,
                                            _apj->Jsi_a_ptr(), _apj->Jsa_a_ptr(),
                                            _apj->Jdia_aa_ptr(), _apj->Jdai_aa_ptr(),
                                            _apj->Jdij_aa_ptr(), _apj->Jdab_aa_ptr(),
                                            weight, mei_adj, tei_adj, local_e_aaaa);
        } else {
          throw formic::Exception("fomric::RGemPow::local_orb_opt_derivs currently only implemented for jagp");
        }

        // compute bb/bbbb contribution to energy and derivatives
        // note that this is done by exchanging alpha and beta and using the alpha code
        if (_apj) {
          this->loc_en_integral_der_aa_aaaa(_nap, _norb-_nap, bo2s, ao2s, bu2s, au2s,
                                            current_im.t(), current_pm.t(),
                                            _apj->Jsi_b_ptr(), _apj->Jsa_b_ptr(),
                                            _apj->Jdia_bb_ptr(), _apj->Jdai_bb_ptr(),
                                            _apj->Jdij_bb_ptr(), _apj->Jdab_bb_ptr(),
                                            weight, mei_adj, tei_adj, local_e_bbbb);
        } else {
          throw formic::Exception("fomric::RGemPow::local_orb_opt_derivs currently only implemented for jagp");
        }

        // compute aabb contribution to energy and derivatives
        if (_apj) {
          this->loc_en_integral_der_aabb(_nap, _norb-_nap, ao2s, bo2s, au2s, bu2s,
                                         current_im, current_pm,
                                         _apj->Jsi_a_ptr(),
                                         _apj->Jsa_a_ptr(),
                                         _apj->Jsi_b_ptr(),
                                         _apj->Jsa_b_ptr(),
                                         _apj->Jdia_aa_ptr(),
                                         _apj->Jdij_ab_ptr(),
                                         _apj->Jdia_ab_ptr(),
                                         _apj->Jdai_ab_ptr(),
                                         _apj->Jdab_ab_ptr(),
                                         _apj->Jdia_bb_ptr(),
                                         weight, tei_adj, local_e_aabb);
        } else {
          throw formic::Exception("fomric::RGemPow::local_orb_opt_derivs currently only implemented for jagp");
        }

        // compute the non-excitation component of the derivatives
        this->non_excite_en_integral_der_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s, formic::unity(S()) * weight, mei_adj, tei_adj);
        this->non_excite_en_integral_der_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s, formic::unity(S()) * weight, mei_adj, tei_adj);

        // combine the components of the local energy
        local_e =   local_e_aaaa
                  + local_e_bbbb
                  + local_e_aabb
                  + this->non_excite_energy_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s)
                  + this->non_excite_energy_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s)
                  + _ham.e_cor_nuc() * formic::unity(S());

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Returns the Hamiltonian contraction with a derivative ratio vector.  Note that the first
      ///        element of the vector corresponds to the no-derivative term.
      ///
      /// \param[in,out]  dfc      the derivative function coefficients.  size == 1+this->n_der()
      ///
      /// \return  sum_i <n|H|Psi_i> dfc_i / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S der_hc(const S * const dfc) {

        // if requested, use the adjoint system
        if (_der_hc_by_adj)
          return this->der_hc_autodiff(dfc);

        formic::start_timer("der_hc");

        // compute local energy
        const S loc_e = this->local_energy();
        //const S loc_e_p = this->partial_local_energy();

        // compute the contribution from the term that does not involve a derivative
        S retval = dfc[0] * loc_e;

        // pairing matrix derivatives
        if (!_fixed_pm) {

          // get a pointer to the pairing matrix derivative coefficients
          const S * const pm_dfc = dfc + 1 + this->pm_pos();

          // compute the contraction of the inverse matrix with the derivative coefficients
          formic::fqmc::rgp_build_occupied_matrix(_nap,
                                                  _norb,
                                                  pm_dfc,
                                                  &_config.a_occ_to_site().at(0),
                                                  &_config.b_occ_to_site().at(0),
                                                  &_oic.at(0));
          formic::xgemm('N', 'N', _nap, _nap, _nap, -formic::unity(S()), &_oic.at(0), _nap, &_oim.at(0),   _nap, formic::zero(S()), &_swork.at(0), _nap);
          formic::xgemm('N', 'N', _nap, _nap, _nap,  formic::unity(S()), &_oim.at(0), _nap, &_swork.at(0), _nap, formic::zero(S()), &_oic.at(0),   _nap);

          // compute contraction of derivative coefficients with derivatives of the regular determinant
          {
            S dc_reg_det = formic::zero(S());
            for ( int i = 0; i < _nap; i++)
              for ( int j = 0; j < _nap; j++)
                dc_reg_det += _oim[_nap*j+i] * pm_dfc[ _config.a_occ_to_site()[i] * _norb + _config.b_occ_to_site()[j] ];
            retval += dc_reg_det * loc_e;
          }

          // prepare derivative contracted new row, new column, and unoccupied pairing matrices
          this->prep_row_col_unocc_mat(pm_dfc, _Zia, _Yjb, _Qab);

          const int no = _nap;
          const int nv = _norb - _nap;

          // prepare derivative contracted inverse matrix and its transpose
          for (int j = 0; j < no; j++)
          for (int i = 0; i < no; i++) {
            _Ckl[i+no*j] = _oic[i*no+j];
            _Clk[j+no*i] = _oic[i*no+j];
          }

          // add aa and aaaa contraction terms
          formic::fqmc::rgp_aaaa_der_hc(no, nv, &_Tia_aa[0], &_TVBjb_aaaa[0], &_TVBaj[0], &_Bkl[0], &_Ckl[0],
                                        &_Via[0], &_Zia[0], &_ZBia[0], &_VCia[0], &retval);

          // add bb and bbbb contraction terms
          formic::fqmc::rgp_aaaa_der_hc(no, nv, &_Tia_bb[0], &_TVBjb_bbbb[0], &_TUBaj[0], &_Blk[0], &_Clk[0],
                                        &_Ujb[0], &_Yjb[0], &_BYjb[0], &_CUjb[0], &retval);

          // add aabb contraction terms
          formic::fqmc::rgp_aabb_der_hc(no, nv, &_Tiajb_aabb[0], &_Via[0], &_Zia[0], &_Ujb[0], &_Yjb[0],
                                        &_Pab[0], &_Qab[0], &_Ckl[0], &_Clk[0],
                                        &_VBia[0], &_VCia[0], &_ZBia[0], &_BUjb[0], &_CUjb[0], &_BYjb[0],
                                        &_VBUab[0], &_TBab[0], &_TVBjb_aabb[0],
                                        &_dVBUab[0], &_TCab[0], &_TBUia[0], &retval);

        }

        // jastrow factor derivatives
        if (_apj) {

          // get number of ocupied and unoccuied sites
          const int no = _nap;
          const int nv = _norb - _nap;

          // get references to occupation arrays
          const std::vector<int> & ao2s = _config.a_occ_to_site();
          const std::vector<int> & bo2s = _config.b_occ_to_site();
          const std::vector<int> & au2s = _config.a_unocc_to_site();
          const std::vector<int> & bu2s = _config.b_unocc_to_site();

          // initialize jastrow derivative intermediates
          _apj->init_der_intermediates(dfc + this->jf_pos());

          // transform the aa and aaaa integrals by the Jastrow factor derivatives
          _apj->der_trans_ints_aa_aaaa(&_Tia_aa.at(0), &_Tiajb_aaaa.at(0));

          //// !!! removing aa oei and aaaa tei integrals as they are now being handled through the adjoint system
          //std::fill(_Tia_aa.begin(), _Tia_aa.end(), formic::zero(S()));
          //std::fill(_Tiajb_aaaa.begin(), _Tiajb_aaaa.end(), formic::zero(S()));

          // compute aa and aaaa one and two electron excitation contribution
          formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_aa[0], &_Tiajb_aaaa[0], &_Bkl[0], &_Via[0],
                                        &_VBia[0], &_TVBjb_aaaa[0], &_TVBaj[0], &retval);

          // transform the bb and bbbb integrals by the Jastrow factor derivatives
          _apj->der_trans_ints_bb_bbbb(&_Tia_bb.at(0), &_Tiajb_bbbb.at(0));

          //// !!! removing bb oei and bbbb tei integrals as they are now being handled through the adjoint system
          //std::fill(_Tia_bb.begin(), _Tia_bb.end(), formic::zero(S()));
          //std::fill(_Tiajb_bbbb.begin(), _Tiajb_bbbb.end(), formic::zero(S()));

          // compute bb and bbbb one and two electron excitation contribution
          formic::fqmc::rgp_aaaa_energy(no, nv, &_Tia_bb[0], &_Tiajb_bbbb[0], &_Blk[0], &_Ujb[0],
                                        &_BUjb[0], &_TVBjb_bbbb[0], &_TUBaj[0], &retval);

          // transform the aabb integrals by the Jastrow factor derivatives
          _apj->der_trans_ints_aabb(&_Tiajb_aabb.at(0));

          //// !!! removing aabb tei integrals as they are now being handled through the adjoint system
          //std::fill(_Tiajb_aabb.begin(), _Tiajb_aabb.end(), formic::zero(S()));

          // compute aabb contribution
          formic::fqmc::rgp_aabb_energy(no, nv, &_Tiajb_aabb[0], &_Via[0], &_Ujb[0], &_Pab[0],
                                        &_Bkl[0], &_Blk[0], &_VBia[0], &_BUjb[0],
                                        &_VBUab[0], &_TBab[0], &_TVBjb_aabb[0], &retval);

          // add non-excitation terms
          retval += _apj->Zne() * this->non_excite_energy_aa_aaaa_aabb(ao2s, bo2s, au2s, bu2s);
          retval += _apj->Zne() * this->non_excite_energy_aa_aaaa_aabb(bo2s, ao2s, bu2s, au2s);

          // add nuclear repulsion term
          retval += _apj->Zne() * _ham.e_cor_nuc();

        }

        formic::stop_timer("der_hc");

        // return the overall contraction
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Returns the identity contraction with a derivative ratio vector.  Note that the first
      ///        element of the vector corresponds to the no-derivative term.
      ///
      /// \param[in,out]  dfc      the derivative function coefficients.  size == 1+this->n_der()
      ///
      /// \return  sum_i <n|I|Psi_i> dfc_i / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S der_ic(const S * const dfc) {

        // non-derivative term
        S retval = dfc[0];

        // pairing matrix terms
        if ( !_fixed_pm ) {
          const S * const der_coeff = dfc + 1 + this->pm_pos();
          retval += formic::fqmc::rgp_det_ratio_deriv(_nap,
                                                      _norb,
                                                      0,
                                                      (int *)0,
                                                      (int *)0,
                                                      _pm,
                                                      &_oim.at(0),
                                                      der_coeff,
                                                      der_coeff,
                                                      &_config.a_occ_to_site().at(0),
                                                      &_config.b_occ_to_site().at(0),
                                                      &_config.a_site_to_occ().at(0),
                                                      &_config.b_site_to_occ().at(0),
                                                      &_swork.at(0));
        }

        // jastrow factor terms
        if (_apj)
          retval += _apj->der_ic(dfc + this->jf_pos()) - dfc[this->jf_pos()];

        // return the overall contraction
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief copies the coefficient's variables to the supplied vector
      ///
      /// \param[out]     vars     on exit, the variables. size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void to_vector(S * const vars) {

        if (!_fixed_pm)
          formic::xcopy(_norb*_norb, _pm, 1, vars + this->pm_pos(), 1);

        if (_apj)
          _apj->to_vector(vars + this->jf_pos());

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief copies the coefficient's variables from the supplied vector
      ///
      /// \param[in]      vars     the variables. size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void from_vector(const S * const vars) {

        if (!_fixed_pm)
          formic::xcopy(_norb*_norb, vars + this->pm_pos(), 1, _pm, 1);

        if (_apj)
          _apj->from_vector(vars + this->jf_pos());

        // clear the exact sampling vectors
        this->clear_exact_sample_vecs();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief reads in all the coefficient object's variables (even the "fixed" ones)
      ///
      /// \param[in]      vars     the variables. size == this->n_var()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void all_from_vector(const S * const vars) {
        
        // get pairing matrix
        formic::xcopy(_norb*_norb, vars + this->pm_pos(false), 1, _pm, 1);

        // get jastrow variables
        if (_apj)
          _apj->all_from_vector(vars + this->jf_pos(false));

        // clear the exact sampling vectors
        this->clear_exact_sample_vecs();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns a vectorization of all the coefficient's variables (even the "fixed" ones)
      //
      /// \param[out]     vars     vector in which to store the variables. size == this->n_var()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void all_to_vector(S * const vars) const {

        // get pairing matrix
        formic::xcopy(_norb*_norb, _pm, 1, vars + this->pm_pos(false), 1);

        // get jastrow variables
        if (_apj)
          _apj->all_to_vector(vars + this->jf_pos(false));

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  updates the coefficient object with a rotation of the orbital basis
      ///
      /// \param[in]      norb     the number of orbitals
      /// \param[in]      orot     the orbital rotation matrix (new/old basis labelled by rows/cols).
      ///                          size == norb*norb
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void apply_orbital_rotation(const int norb, const S * const orot) {
        throw formic::Exception("still need to implement rgp's apply_orbital_rotation function");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief prints the coefficient's variables to the output file
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void print_variables() {

        if ( formic::mpi::rank() == 0 ) {
          formic::of << "printing restricted geminal power pairing matrix:" << std::endl;
          for (int i = 0; i < _norb; i++) {
            for (int j = 0; j < _norb; j++)
              formic::of << " " << formic::format_number("%14.8f", _pm[_norb*i+j]);
            formic::of << std::endl;
          }
          formic::of << std::endl;
        }

        if (_apj)
          _apj->print_variables();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of nonzero matrix elements in the sparse preconditioner
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      int n_pc_nonzero() {

        // get number of rows corresponding to one-orbital variables
        int nii = 0;
        if (!_fixed_pm)
          nii += _norb;
        if (_apj)
          if (!_apj->ab_fixed())
            nii += 4 * _norb;

        // get number of rows corresponding to two-orbital variables
        int nij = 0;
        if (!_fixed_pm)
          nij += _norb*_norb - _norb;
        if (_apj) {
          if (!_apj->aa_fixed())
            nij += _apj->naa_var();
          if (!_apj->bb_fixed())
            nij += _apj->naa_var();
          if (!_apj->ab_fixed())
            nij += 4 * (_norb*_norb - _norb);
        }

        // compute number of nonzero preconditioner elements
        int retval = 1 + 2*this->n_der(); // count the full first row and column
        if (!_fixed_pm)
          retval += nii * (2*_norb-1) + nij * (4*_norb-4);
        if (_apj) {
          if (!_apj->aa_fixed())
            retval += nii * 4*(_norb-1) + nij * 4*(2*_norb-3);
          if (!_apj->bb_fixed())
            retval += nii * 4*(_norb-1) + nij * 4*(2*_norb-3);
          if (!_apj->ab_fixed())
            retval += nii * 4*(2*_norb-1) + nij * 4*(4*_norb-4);
        }

        // return the total
        return retval;

      }

    private:

      void pc_row_inds_ii(const int i_in, int * & n_ptr, int * & ind_ptr) {

        struct local_funcs {

          static void jf_aa_ii(int & base_ind,
                               const int naa_var,
                               int & n,
                               const int norb,
                               const int i_in,
                               int * & ind_ptr) {

            int tot = 0;
            for (int j = 1; j < i_in; j++)
              tot += j;

            for (int i = 0; i < i_in; i++) {
              const int jf_ind = base_ind + 4*(tot+i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }

            for (int j = i_in+1; j < norb; j++) {
              tot += j-1;
              const int jf_ind = base_ind + 4*(tot+i_in);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }

            base_ind += naa_var;
            n += 4*(norb-1);

          }

        };

        // zeroeth column is always touched
        int n = 1;
        *ind_ptr++ = 0;

        // pairing matrix columns
        if (!_fixed_pm) {
          n += 2*_norb - 1;
          const int base_ind = 1 + this->pm_pos();
          for (int i = 0; i < i_in; i++)
            *ind_ptr++ = base_ind + _norb*i + i_in;
          for (int j = 0; j < _norb; j++)
            *ind_ptr++ = base_ind + _norb*i_in + j;
          for (int i = i_in + 1; i < _norb; i++)
            *ind_ptr++ = base_ind + _norb*i + i_in;
        }

        if (_apj) {

          int base_ind = 1 + this->jf_pos();

          // aa jastrow columns
          if (!_apj->aa_fixed())
            local_funcs::jf_aa_ii(base_ind, _apj->naa_var(), n, _norb, i_in, ind_ptr);

          // bb jastrow columns
          if (!_apj->bb_fixed())
            local_funcs::jf_aa_ii(base_ind, _apj->naa_var(), n, _norb, i_in, ind_ptr);

          // ab jastrow columns
          if (!_apj->ab_fixed()) {
            for (int j = 0; j < i_in; j++) {
              const int jf_ind = base_ind + 4*(_norb*j + i_in);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            for (int i = 0; i < _norb; i++) {
              const int jf_ind = base_ind + 4*(_norb*i_in + i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            for (int j = i_in + 1; j < _norb; j++) {
              const int jf_ind = base_ind + 4*(_norb*j + i_in);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            base_ind += _apj->naa_var();
            n += 4*(2*_norb-1);
          }

        }

        // record the total number of nozero columns for the row
        // and increment the row pointer
        *n_ptr++ = n;

      }

      void pc_row_inds_ij(int a, int b, int * & n_ptr, int * & ind_ptr) {

        struct local_funcs {

          static void jf_aa_ij(int & base_ind,
                               const int naa_var,
                               int & n,
                               const int norb,
                               const int a,
                               const int b,
                               int * & ind_ptr) {

            int tot = 0;
            for (int j = 1; j < a; j++)
              tot += j;

            for (int i = 0; i < a; i++) {
              const int jf_ind = base_ind + 4*(tot+i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            tot += a;

            for (int j = a+1; j < b; j++) {
              const int jf_ind = base_ind + 4*(tot+a);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
              tot += j;
            }

            // b-1 up to here

            for (int i = 0; i < b; i++) {
              const int jf_ind = base_ind + 4*(tot+i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            tot += b;

            // 2*b-1 up to here

            // the following for loop gives an additional 2*(norb-b-1)
            for (int j = b+1; j < norb; j++) {
              const int jfa = base_ind + 4*(tot+a);
              *ind_ptr++ = jfa + 0;
              *ind_ptr++ = jfa + 1;
              *ind_ptr++ = jfa + 2;
              *ind_ptr++ = jfa + 3;
              const int jfb = base_ind + 4*(tot+b);
              *ind_ptr++ = jfb + 0;
              *ind_ptr++ = jfb + 1;
              *ind_ptr++ = jfb + 2;
              *ind_ptr++ = jfb + 3;
              tot += j;
            }

            // total of 2*norb - 3
            n += 4*(2*norb-3);
            base_ind += naa_var;

          }

        };

        // ensure a is less than b
        if ( a > b )
          std::swap(a,b);

        // ensure a and b are distinct
        if ( a == b )
          throw formic::Exception("input indices must be distinct in fomric::RGemPow::pc_row_inds_ij");

        // zeroeth column is always touched
        int n = 1;
        *ind_ptr++ = 0;

        // pairing matrix columns
        if (!_fixed_pm) {
          n += 4*_norb - 4;
          const int base_ind = 1 + this->pm_pos();
          for (int i = 0; i < a; i++) {
            *ind_ptr++ = base_ind + _norb*i + a;
            *ind_ptr++ = base_ind + _norb*i + b;
          }
          for (int j = 0; j < _norb; j++)
            *ind_ptr++ = base_ind + _norb*a + j;
          for (int i = a+1; i < b; i++) {
            *ind_ptr++ = base_ind + _norb*i + a;
            *ind_ptr++ = base_ind + _norb*i + b;
          }
          for (int j = 0; j < _norb; j++)
            *ind_ptr++ = base_ind + _norb*b + j;
          for (int i = b+1; i < _norb; i++) {
            *ind_ptr++ = base_ind + _norb*i + a;
            *ind_ptr++ = base_ind + _norb*i + b;
          }
        }

        if (_apj) {

          int base_ind = 1 + this->jf_pos();

          // aa jastrow columns
          if (!_apj->aa_fixed())
            local_funcs::jf_aa_ij(base_ind, _apj->naa_var(), n, _norb, a, b, ind_ptr);

          // bb jastrow columns
          if (!_apj->bb_fixed())
            local_funcs::jf_aa_ij(base_ind, _apj->naa_var(), n, _norb, a, b, ind_ptr);

          // ab jastrow columns
          if (!_apj->ab_fixed()) {
            for (int j = 0; j < a; j++) {
              const int jfa = base_ind + 4*(_norb*j + a);
              *ind_ptr++ = jfa + 0;
              *ind_ptr++ = jfa + 1;
              *ind_ptr++ = jfa + 2;
              *ind_ptr++ = jfa + 3;
              const int jfb = base_ind + 4*(_norb*j + b);
              *ind_ptr++ = jfb + 0;
              *ind_ptr++ = jfb + 1;
              *ind_ptr++ = jfb + 2;
              *ind_ptr++ = jfb + 3;
            }
            for (int i = 0; i < _norb; i++) {
              const int jf_ind = base_ind + 4*(_norb*a + i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            for (int j = a+1; j < b; j++) {
              const int jfa = base_ind + 4*(_norb*j + a);
              *ind_ptr++ = jfa + 0;
              *ind_ptr++ = jfa + 1;
              *ind_ptr++ = jfa + 2;
              *ind_ptr++ = jfa + 3;
              const int jfb = base_ind + 4*(_norb*j + b);
              *ind_ptr++ = jfb + 0;
              *ind_ptr++ = jfb + 1;
              *ind_ptr++ = jfb + 2;
              *ind_ptr++ = jfb + 3;
            }
            for (int i = 0; i < _norb; i++) {
              const int jf_ind = base_ind + 4*(_norb*b + i);
              *ind_ptr++ = jf_ind + 0;
              *ind_ptr++ = jf_ind + 1;
              *ind_ptr++ = jf_ind + 2;
              *ind_ptr++ = jf_ind + 3;
            }
            for (int j = b+1; j < _norb; j++) {
              const int jfa = base_ind + 4*(_norb*j + a);
              *ind_ptr++ = jfa + 0;
              *ind_ptr++ = jfa + 1;
              *ind_ptr++ = jfa + 2;
              *ind_ptr++ = jfa + 3;
              const int jfb = base_ind + 4*(_norb*j + b);
              *ind_ptr++ = jfb + 0;
              *ind_ptr++ = jfb + 1;
              *ind_ptr++ = jfb + 2;
              *ind_ptr++ = jfb + 3;
            }
            base_ind += _apj->naa_var();
            n += 4*(4*_norb-4);
          }

        }

        // record the total number of nozero columns for the row
        // and increment the row pointer
        *n_ptr++ = n;

      }

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief computes the index arrays for the sparse preconditioner
      ///
      /// \param[in,out]  pc_n_inds  size 1+this->n_der().
      ///                            On exit, the number of nonzero matrix elements for each row.
      /// \param[in,out]  pc_inds    size this->n_pc_nonzero().
      ///                            On exit, the column indices of each rows' nonzero elements.
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void set_pc_inds(int * const pc_n_inds, int * const pc_inds) {

        int * n_ptr = pc_n_inds;
        int * ind_ptr = pc_inds;

        // first row is full
        *n_ptr = 1+this->n_der();
        for (int i = 0; i < *n_ptr; i++)
          *ind_ptr++ = i;
        n_ptr++;

        // pairing matrix variables
        if (!_fixed_pm) {
          for (int i = 0; i < _norb; i++)
          for (int j = 0; j < _norb; j++)
            if (i == j)
              this->pc_row_inds_ii(i, n_ptr, ind_ptr);
            else
              this->pc_row_inds_ij(i, j, n_ptr, ind_ptr);
        }

        // all pair jastrow
        if (_apj) {

          // aa jastrows
          if (!_apj->aa_fixed())
            for (int j = 0; j < _norb; j++)
            for (int i = 0; i < j; i++)
            for (int k = 0; k < 4; k++)
              this->pc_row_inds_ij(i, j, n_ptr, ind_ptr);

          // bb jastrows
          if (!_apj->bb_fixed())
            for (int j = 0; j < _norb; j++)
            for (int i = 0; i < j; i++)
            for (int k = 0; k < 4; k++)
              this->pc_row_inds_ij(i, j, n_ptr, ind_ptr);

          // ab jastrows
          if (!_apj->ab_fixed())
            for (int j = 0; j < _norb; j++)
            for (int i = 0; i < _norb; i++)
            for (int k = 0; k < 4; k++)
              if (i == j)
                this->pc_row_inds_ii(i, n_ptr, ind_ptr);
              else
                this->pc_row_inds_ij(j, i, n_ptr, ind_ptr);

        }

        if ( n_ptr != pc_n_inds + 1 + this->n_der() )
          throw formic::Exception("bad iteration over n_ptr");

        if ( ind_ptr != pc_inds + this->n_pc_nonzero() )
          throw formic::Exception("bad iteration over ind_ptr");

        int problem = 0;
        for (int i = 0, stop = this->n_pc_nonzero(), max_ind = 1 + this->n_der(); i < stop; i++) {
          problem += ( pc_inds[i] >= max_ind ? 1 : 0 );
          problem += ( pc_inds[i] < 0 ? 1 : 0 );
        }
        if (problem != 0)
          throw formic::Exception("bad value in pc_inds");

        for (int i = 0, max_ind = 1 + this->n_der(); i < max_ind; i++) {
          problem += ( pc_n_inds[i] > max_ind ? 1 : 0 );
          problem += ( pc_n_inds[i] < 1 ? 1 : 0 );
        }
        if (problem != 0)
          throw formic::Exception("bad value in pc_n_inds");

        if ( formic::mpi::rank() == 0 )
          formic::of << "preconditioner indices initialized" << std::endl << std::endl;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief clones the object
      ///
      /// \return  a shared pointer to a new clone of the object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      boost::shared_ptr<formic::CoeffObj<S> > clone() const {

        boost::shared_ptr<formic::CoeffObj<S> > retval( new RGemPow<S>(*this) );
        return retval;

      }

  };

}

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "RGemPow" command, which adds an AGP coefficient object
///          to the wavefunction
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////

static void formic_RGemPow_command_function(formic::Archive & arch, formic::AllData & data) {

  typedef std::complex<double> scdoub;

  // read in and add the restricted geminal power coefficient object to the wavefunction
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    boost::shared_ptr<formic::CoeffObj<double> > c( new formic::fqmc::RGemPow<double>(data.userinp, data.hd, data.wd, arch, false) );
    data.wd.add_coeff_obj(c);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    boost::shared_ptr<formic::CoeffObj<scdoub> > c( new formic::fqmc::RGemPow<scdoub>(data.userinp, data.hc, data.wc, arch, false) );
    data.wc.add_coeff_obj(c);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing RGemPow command.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "JRGemPow" command, which adds an all-pair jastrow AGP
///          coefficient object to the wavefunction
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////

static void formic_JRGemPow_command_function(formic::Archive & arch, formic::AllData & data) {

  typedef std::complex<double> scdoub;

  // read in and add the restricted geminal power coefficient object to the wavefunction
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    boost::shared_ptr<formic::CoeffObj<double> > c( new formic::fqmc::RGemPow<double>(data.userinp, data.hd, data.wd, arch, true) );
    data.wd.add_coeff_obj(c);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    boost::shared_ptr<formic::CoeffObj<scdoub> > c( new formic::fqmc::RGemPow<scdoub>(data.userinp, data.hc, data.wc, arch, true) );
    data.wc.add_coeff_obj(c);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing JRGemPow command.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Add restricted geminal power related commands to the formic command system
///
///////////////////////////////////////////////////////////////////////////////////////////////////

void formic::init_rest_gem_pow_commands() {

  // add command to construct an AGP coefficient object
  formic::add_command("RGemPow", &formic_RGemPow_command_function);

  // add command to construct an all-pair jastrow AGP coefficient object
  formic::add_command("JRGemPow", &formic_JRGemPow_command_function);

}
